Dynamics modeling method dedicated to system models with open- and closed-loop structures subjected to kinematic and task based constraints

The paper presents a development of an automated computational procedure for constrained dynamics (CoPCoD) dedicated to derivation of dynamical models of mechanical systems, e.g. manipulators, both ground and space or mobile robotic systems, which can be composed of rigid and flexible links. They may be subjected to constraints, which are referred to as programmed, which may come from performance requirements, e.g. work or services a system is dedicated to. The CoPCoD structure offers systematic modeling of either open or closed loop constrained structures and results in computationally efficient numerical dynamical equations derivation. The CoPCoD development has its background in the generalized programmed motion equations (GPME) algorithm developed successfully for rigid models of mechanical systems subjected to high order nonholonomic constraints, however, the GPME was not fully automated for computer equation derivations. The two main motivations that underlie the CoPCoD development are to extend the GPME to flexible system models and make it fully automated to computer equation derivations for large classes of systems. The paper presents the general scheme of the CoPCoD architecture and examples which illustrate its application advantages.

Abstract The paper presents a development of an automated computational procedure for constrained dynamics (CoPCoD) dedicated to derivation of dynamical models of mechanical systems, e.g. manipulators, both ground and space or mobile robotic systems, which can be composed of rigid and flexible links. They may be subjected to constraints, which are referred to as programmed, which may come from performance requirements, e.g. work or services a system is dedicated to. The CoPCoD structure offers systematic modeling of either open or closed loop constrained structures and results in computationally efficient numerical dynamical equations derivation. The CoPCoD development has its background in the generalized programmed motion equations (GPME) algorithm developed successfully for rigid models of mechanical systems subjected to high order nonholonomic constraints, however, the GPME was not fully automated for computer equation derivations. The two main motivations that underlie the CoPCoD development are to extend the GPME to flexible system models and make it fully automated to computer equation derivations for large classes of systems. The paper presents the general scheme of the CoPCoD architecture and examples which illustrate its application advantages. Stiffness, damping and viscous friction coefficients U Programmed and kinematic constraint equations U k ; U p Kinematic and programmed constraints equations rfe l Rigid fininte element related to a flexible link sde s Spring-damping element related to the support sde l Spring-damping element related to a flexible link 1 Introduction

Background and literature survey
The literature on various aspects of research on rigid and flexible multibody systems is very rich. There is a lot of papers about writing, solving and simulation of dynamic equations of rigid and flexible multibody systems; see e.g. [1,2] and references there. Also, there are multiple papers on designing control techniques for such systems. Depending upon the purpose of modeling, we can distinguish methods aimed to build validation models including nonlinear terms [3,4], methods to derive preliminary design models, which due to their simplicity are attractive to control [5]. The control dedicated models often arise as various versions of linearized equations, small displacement assumptions are adopted for flexible multibody systems and link flexibility is approximated by discretized models using finite element methods or assumed mode methods; see e.g. [6]. The dynamics of the system is usually derived using the Euler-Lagrange formulation. As indicated in [7], there are cases, e.g. in space manipulation during debris capture, when there is a need to have a model for each sub-structure of a system such that it is valid for arbitrary boundary conditions. Among models that take advantage of the substructure approaches, the transfer matrix method gained popularity [8]. Later, this approach was merged with the finite-element method and presented as the finite-element transfer matrix method [9]. Some extensions of this method were also developed for control design purposes and for rigid and flexible multibody systems [10]. However, the transfer matrix method has some limitations due to the inversion of sub-matrices, which are not always square or invertible, depending on the boundary conditions. It cannot be directly applied to tree-like structures, like a spacecraft composed of a rigid base and flexible appendages on it.
In many works various methods for derivation of linear or linearized models of flexible multibody system dynamics are presented. For example in [7], the approach is based on the two-port model of each body allowing the whole system model to be built by connecting the inputs/outputs of each body model. Boundary conditions of each body are taken into account through inversion of some input-output channels of its two-port model. This one model of the link can be built into an open or closed-loop mechanism. The method presented in [7] fills the gap between the transfer matrix method and the effective mass-inertia method.
The significant group of dynamic models in engineering applications origins from the multibody dynamics and control of rigid and flexible systems for which Newton-Euler and Lagrange based dynamics derivation methods are used. Among them, there are special methods dedicated to constrained systems in both analytical and computational approaches. As reported in [11], constrained system dynamical models constitute a significant class of models essential to engineering applications. They are used so videly, e.g. in motion analysis in robotics for multi-link industrial manipulators and robots, for ground, underwater, aerospace and space vehicle system dynamics, control, and performance analysis for variety of their missions, that computational methods for derivation of their dynamical models, and their solution methods, became separate research areas, see e.g. works [12,13] and references there. A lot of application examples for multibody system dynamics models can be found in [14]. Reviewing derivation methods of dynamic models for mechanical systems, the first observation is that they are for systems subjected to position and first order constraints, which are material ones; see e.g. [15,16] and references there. The task based constraints, control and performance constraints and other requirements available in equation forms are not merged into these models. Consequently, these models are derived using classical analytical dynamics methods, i.e. based upon the Lagrange approach and its modifications. Another approach to dynamics modeling and control approach is based upon the framework developed by Udwadia and Kalaba [17,18]. This is an elegant apparatus for deriving constrained dynamics based upon the classical mechanics approach and it serves systems subjected to position and kinematic constraints. Also, there are many specialized computational packages for generation of constrained system dynamical models. Generally, they are either Newton-Euler equations or Lagrange equations based and dedicated to rigid multibody systems; see e.g. [19,20]. The consequence of the derivation method applied to obtain constrained system models is that the resulting equations of motion may not be convenient for their final applications, e.g. they contain Lagrange multipliers which need to be eliminated for most model based controller designs. The second observation is that some of the proposed methods of motion equations derivation are specific and may be challenging for applications for systems with flexible components; see e.g. [21]. The good illustrations of the two observations mentioned above are works [22][23][24]. These works have had significant impact upon the insight into the constrained system structures and for the subsequent control design, like [24] where geometric properties of dynamic models were shown to be the key to model properties from stability and control points of view. However, constraint kinds considered there and presented results were for material constraints, position and kinematic ones. Also, new challenges in dynamic modeling and control designs emerged from interdisciplinary or specifically dedicated models as well as from requirements for system performance or their service mission demands. The new challenges were brought, for example, by underwater vehicles and spacecraft, see e.g. [25]. For the latter ones, long and lightweight structures require link flexibility incorporation into dynamic models as well as mission service and control requires formulation of task based constraints and some of them can be easily formulated by constraint equations. Usually such constraints, if formulated, are at position levels; see e.g. [26,27]. Also, one more challenge arises for spacecraft modeling and simulation; this is to keep conserved quantities of motion like the linear and angular momentum. The linear and angular momenta can be considered the kinematic constraints and incorporated into constrained dynamics [28].

Formulation of the problem of interest for this investigation
In most dynamics derivation methods, of either approaches, systems subjected to position and velocity constraints are considered and the Lagrange approach is used. Other constraints are not handled but there are the ones, which are generated by designers, control engineers or system operators that have to be executed during a system operation, see e.g. [29][30][31][32]. For systems subjected to such constraints the preliminary results for obtaining computationally efficient dynamical models were presented in [11,33]. Such constraints at the position level, referred to as programmed, were incorporated into a rigid-flexible system dynamics for open loop structures. The analytical dynamics background for the constrained dynamic models are works [34,35], where the generalized programmed motion equations (GPME) are derived. The GPME proved to be effective to derive constrained dynamics, referred to as reference dynamics, for rigid system models in various parameter settings like generalized or quasi-coordinates. The GPME were employed to design a control strategy architecture for tracking predefined motions. Also, the GPME approach can be applied to variable mass and configuration systems [36]. However, The GPME depend heavily upon analytical dynamics calculations and the merit of this paper is to extend to multi-degree of freedom systems whose dynamics require computation.

Motivation and contribution of this study
This paper is motivated by two factors, being quite critical in designing constrained dynamics and controllers for programmed constraints execution. The first one is the specification of programmed constraints as demands for some work regimes, e.g. velocity of a manipulator end effector has to be modulated during the manipulator work. From the control points of view, incorporation of these constraints into the system dynamics can be advantageous. The second motivation refers to system structures, e.g. rigid or flexible parts with open-and closed-loop structures, for which one modeling method would be welcome. Finally, dynamics modeling should be proactive not reactive what means that its destination should be clear before but not after modelling [37]. The main contribution of the paper is development of the automated computational procedure for constrained dynamics (CoPCoD) of the flexibly supported systems. In contrast to the authors' earlier work, programmed constraints are defined as rheonomic. The presented approach is based on the GPME and preserves its properties, however is fully computationally automated. The basic one is that the constrained dynamics it produces is in the reduced state form, i.e. constraint reaction force are eliminated at the derivation level but not afterwards as in the Lagrange approach. This is the essential advantage of the CoPCoD which may serve directly to constrained motion analysis as well as to control design. The CoPCoD offers automation of the constrained dynamics generation for rigid and flexible system models with open-and closed-loop structures, and possibility of incorporation into these dynamics any position or first order material or non-material constraints, also these related to operation regimes. The credibility of the CoPCoD method is built by its development on the GPME's. The GPME proved to be effective to constrained rigid multibody systems dyanmics and control. The CoPCoD enjoys all the GPME's properties and is substantially extended to flexible structures with open and closed loops.
The paper presents the general structure of the CoPCoD and demonstrates the ease of its application in cases of complex kinematic structures, e.g. a tree structure with open and closed loop chains. Examples illustrate how the CoPCoD structure supports modeling and computation of complex structures.

Organization of the paper
The paper is organized as follows. After presenting the background and motivation of the paper, the general algorithm of CoPCoD for open-and closed-loop serial kinematic chains subjected to programmed constraints is presented in Sect. 2. In this section, models of a flexible base support, a flexible link and friction in joints are described.  spring-damping elements (sde).

Generalized coordinates and homogeneous transformation matrices
Motion of a system can be defined by a vector of generalized coordinates as follows: where: n dof is number of generalized coordinates which describe a system motion, n dof ¼ In the case, when the system reference motion is analyzed, i.e. according to constraints imposed upon it, the vector of generalized coordinates is divided into independent and dependent ones, i.e.
where i i c ; i d c are indices of the independent and dependent coordinates, respectively. Motion of a link a; i ð Þj a2 c;c 1 ;c 2 f g with respect to the global reference frame F ð0Þ is described by the following vector: where: n ðaÞ l is the number of links (rigid and flexible) in chain a; n ðaÞ l ¼ n ðaÞ r l þ n ðaÞ f l ;q ða;iÞ is a vector defining motion of link a; i ð Þ with respect to a preceding link a; i À 1 ð Þand q ða;0Þ 2 ;: When the link a; i ð Þis treated as flexible, its motion with respect to the preceding link is described by a vector: whereq ða;i;rÞ is a vector defining motion of rfe a; i; r ð Þ with respect to the preceding rfe a; i; r À 1 ð Þ :. Homogeneous transformation matrices from a local frame of a link to the global reference frame F ð0Þ can be obtained in the following way: where: T ðbÞ is the transformation matrix from a frame F ðbÞ fixed at the flexibly supported base to the global reference frame F ð0Þ andT ða;jÞ is the transformation matrix from a local frame of the link a; i ð Þ to a local frame of a preceding link in a kinematic chain.
The homogeneous transformation matrices of rfes from their local frames F ða;i;rÞ to the global reference frame F ð0Þ are determined by:

Generalized Programmed Motion
Equations (GPME) for systems subjected toposition and velocity constraints The constrained dynamics, referred to as the reference dynamics, is obtained based upon a dedicated method to constrained dynamics generation, i.e. upon the Generalized Programmed Motion Equations (GPME). The method developed for constrained rigid system structures subjected to arbitrary order nonholonomic constraints, which can specify motion demands and are called programmed, is detailed in [34,35]. It enables verification of the system behavior when it is subjected to programmed constraints, e.g. vibrations that may accompany the desired motions specified by these constraints. It enables then defining and analyzing desired motions for performing servicing tasks. As detailed in [34,35], the constraint equations imposed as control goals on system performance or as service tasks and can be presented in a general form as: B t; q; _ q; :::; q ðpÀ1Þ q ðpÞ þ s t; q; _ q; :::; q ðpÀ1Þ ¼ 0; where: q ðpÞ is a n dof -dimensional state vector, B is n dof ! n d c , and s is a n d c -dimensional vector, n d c is the number of dependent generalised coordinates. The constraints can be material for p ¼ 0; 1, or programmed for p ! 1. The programmed constraints are imposed by a designer or a control engineer to obtain a system desired performance, e.g. they can be imposed upon acceleration p ¼ 2 or jerk p ¼ 3, as well as for desired trajectories with p ¼ 1. The constraints form (7) is the generalized constraint formulation and it encompasses the classical analytical constraint concept.
The GPME for rigid body models subjected to constraints (7) have the form: ð8:1Þ B t; q; _ q; :::; q ðpÀ1Þ q ðpÞ þ s t; q; _ q; :::; q ðpÀ1Þ ¼ 0; ð8:2Þ Þ is a vector of centrifugal forces, g q ð Þ is a vector of gravity forces, Þ is a vector of generalized forces vector resulting from external actions which are not controls. The set of Eqs. (8.1) and (8.2) can be transformed to the form that uses directly functions for the kinetic and potential energies, as well as explicit forms of external forces on a system. The GPME modified for flexible open-loop structures subjected to position constraints are detailed in [33].
For the development of a general structure of the automated CoPCoD, we take advantage of the form of the GPME that uses the energy and force functions. Also, the actual form of the CoPCoD can incorporate position and velocity constraints, both material and non-material into system dynamics. From engineering application points of view such constraints are the most significant ones. Modifications that are dedicated to motion description of rigid and flexible open and closed-loop substructures have to be substantiated. The modified new structure of the GPME for position and velocity constraints is introduced as follows where now the function R 1 that collects the system energies and forces is composed of for CLKC; are kinetic energy terms of the flexibly supported base and the chains, p;g is the potential energy of the base which can be calculated as a sum of spring deformations energy and the potential energy of the gravity forces, p;f l is the potential energy of the chain which can be expressed as a sum of spring deformations energy of flexible links and the potential energy of the gravity forces of links, Thus, in order to derive the CoPCoD, it is required to determine the kinetic and potential energies, the Rayleigh dissipation functions and external forces for components of kinematic chains, as well as formulate the constraint equations.

Kinetic energy of the system
The kinetic energy of the system structure can be presented in the following forms: The time derivatives of the transformation matrices are determined by the following formulas:  The flexible support of the base is modelled by means of spring-damping elements as shown in Fig. 2.
Positions of the supports are described with respect to the base fixed frame. It is assumed that deformations of the spring-damping elements are small and the transformation matrix can be linearized. The spring deformation energy and the Rayleigh dissipation function of the spring-damping elements are determined according to the following formulas: It can be seen that deformations and deformation velocities of sde l ða; i; jÞ are taken directly from the generalized coordinates vector and its derivative.

Generalized forces resulting from actions of friction torques and forces
The LuGre friction model [39] is used to describe friction phenomenon in joints. Planar models of revolute and prismatic joints with friction are presented in Fig. 3. The generalized forces vector of a chain a resulting from including friction in joints can be written as follows:  The joint forces can be determined, in each integration step of the dynamics equations of motion, using the recursive Newton-Euler algorithm.
Friction coefficients are calculated using the LuGre friction model from the following formula:

Kinematic and programmed constraint equations
The GPME (9.1) needs to be supplemented by constraint equations composed of n k kinematic and n p programmed constraint equations which can be written in the following form: where U k q ð Þ; U p q ð Þ are kinematic and programmed constraints equations, g p ðtÞ is a vector of known time functions. The kinematic constraints come from closed chains modelling and are material constraints.
Time derivatives of the constraint Eqs. (23) lead to: ð24:1Þ where: U q;i c ; U q;d c are obtained the from the global constraint matrix U q ¼ U k;q U p;q ! by selecting rows and columns corresponding to indices of the independent and dependent variables in vector q.
where:  ;    ( Equation (26) form a set of n dof differential algebraic-equations. It can be seen that the presented approach allows to obtain the dynamics free of the constraint reaction forces. Also, CoPCoD is structured into block architecture such that different systems can be modelled by selecting options provided by its blocks.

Examples of application of CoPCoD applications
In this section, two sample systems representing OLKC and CLKC structures are analysed. We demonstrate how the CoPCoD algorithm enables performing efficient numerical simulations. They provide information about the system model behaviour when the programmed constraints are imposed upon motion. It enables, in turn, constrained motion planning.

The OLKC manipulator model
As the first example, the programmed motion of the spatial OLKC manipulator model is analyzed. The manipulator possesses 5 degrees of freedom with respect to the base frame {b} (Fig. 4). It is assumed that link (1,2) can be flexible and flexibility of the ground is modelled by means of the n sde s springdamping elements, where n sde s ¼ P n e i¼1 n ðiÞ sde s : The local frames assigned to each link and generalized coordinates of the manipulator are presented in Fig. 5  It is assumed that the end-effector E has to track the rectangle trajectory defined in y ðPÞ z ðPÞ plane (Fig. 4) and it moves between adjacent corners of the rectangle within the scheduled time Dt ðEÞ ¼ 1 s. The change of the end-effector coordinates in this motion is described by a 5 th order polynomial: where: By differentiating (29) with respect to time, one obtains the constraint equations at the velocity and acceleration levels:     When the model with a fixed base is considered, additional constraint equations need to be defined as Finally, equations of motion (26) for the OLKC manipulator model have the form where a; b-the Baumgarte method parameters for dof constraint matrix resulting from (30.1), (30.2), I 6Â6 -identity matrix. In the case when the base support is flexible the constraint (32) and its derivatives are not included into (33).
In Table 1 mass parameters of the manipulator are provided. Friction parameters, stiffness, damping and other parameters are gathered in Table 2.
The influence of the base support stiffness, link's flexibility and friction in joints in the programmed motion performed by the manipulator is analyzed next. Also, forces acting in the joints are inspected. To facilitate the description of subsequent simulation results, the following notation shown in Fig. 6 is introduced. The 4th Runge-Kutta scheme with a constant step size equal to 10 À5 s is applied to integrate the dynamic equations of motion. The Baumgarte method parameters for solution stabilization are selected to be a ¼ 10 3 ; b ¼ 10 2 :  Figure 7 shows the trajectory and time courses of the end-effector. The starting point of the effector motion differs from the assumed point A. The elimination of disturbances resulting from an arbitrary initial configuration of the manipulator end effector is obtained during its transition from the starting point * to B. After this time, the effector moves along the edges of the ABCD rectangle. The effector is stopped at the corner points and then moves to the next point. Figure 8 shows the time course of the magnitude of constraint violations at the displacement level for all analyzed cases. A logarithmic scale on the vertical axis is introduced to increase the resolution of the presented results. Analyzing the results, it can be seen that during the motion from point B to points C, D, and A, violations of the programmed constraints do not exceed the value of 5 Â 10 À3 m. The constraint Fig. 7 The calculated trajectory of the end-effector and time courses of its coordinates Fig. 6 The scheme of notation applied in the simulation results description  Time courses of the joint coordinates are presented in Fig. 9. These results show that the most significant influence on the manipulator motion comes from the flexibility of the base support and friction in joints.
Comparing the plots for the manipulator with the rigid link (1,2) and ideal joints (S r L r J nf ) with those for the flexible link (1,2) and frictionless joints (S r L fl J nf ), one can conclude that they are close to each other. Figure 10 displays plots of the time courses of force magnitudes acting in the joints. It can be seen that In order to analyze the values of the axial and normal forces acting in the joints in the respected cases, the RMS (Root Mean Square) indicator is introduced where: t k is the analysis time, FðtÞ is the time course of the axial or normal forces. Figure 11 allows us to analyze RMS values of the particular components of the joint forces.
Analyzing the RMS values calculated for the axial and normal forces it can be observed that the highest values of the RMS are obtained for the case S r L fl J nf ;, i.e. the link (1,2) is flexible, the support is rigid and friction in joints is neglected. The RMS values determined for this case are up to 8 times bigger than the others. Thus, it can be concluded that the link flexibility, in this case, does not have major impact on the manipulator motion, but may cause a significant increase of the joint loads.
The presented analysis of the system performance in programmed motion enables detailed verification and motion planning for constrained tasks.

The RPSUP mechanism model
The RPSUP (R-revolute, P-prismatic, S-spherical, Uuniversal, P-prismatic joint) mechanism is the second example and it is the CLKC system model (Fig. 12). The goal of the simulation is to track programmed motion by the end effector E, which is a rectangular trajectory located on the plane y ðPÞ z ðPÞ that is inclined to x ð0Þ z ð0Þ plane at the angle a ¼ 45 . This trajectory reflects some work regime on the end effector. The transformation matrix from the frame fPg to the global reference frame has the form: In order to formulate the dynamic equations of motion of the mechanism using the CoPCoD algorithm, first it is necessary to remove the kinematic loop using the cut-joint technique. It is assumed that the cut-joint is located at the revolute joint R connecting the flexible link (1,3) with the slider (2,1). Local coordinate systems are assigned to the links as shown in Fig. 13  Next, applying the CoPCoD algorithm presented in the paper, the dynamic equations of motion of the mechanism are obtained as: Equations (36) are additionally supplemented by closing constraint equations formulated for the revolute cut-joint R. In simulations, as in the case of the manipulator, the influence of the base support flexibility, the coupler flexibility, friction in the joints on the generated programmed motion courses and the forces acting in the joints are analyzed. Mass parameters of the mechanism links applied in simulation tests are presented in Table 3. Other parameters are gathered in Table 2.
Runge-Kutta method is used to integrate dynamic Eqs. (36) with a constant step size h ¼ 10 À5 s. The   applied Baumgarte method coefficient values are a ¼ 500, b ¼ 50. Initial conditions for the dynamic analysis are obtained from solving the statics case. Figures 14 and 15 show the obtained end-effector E trajectory and the corresponding time courses of the end-effector coordinates.
Analyzing the solutions obtained from statics, it can be noticed that the application of the base support flexibility leads to the effector displacement with respect to the point A. Moreover, the coupler flexibility has no major effect on the effector position change at the initial time. Initial arbitrary location of the end effector is eliminated during its transition from the point A to B, and then the mechanism moves according to the programmed constraints. Some distortions of the trajectory for y-component of the end-effector displacement can be observed for the case of S r L r J nf , in which the base support and the coupler are rigid and all joints are treat as ideal. These distortions are even greater at the velocity and acceleration levels (Fig. 16). Figure 17 displays magnitudes of the constraint violations for the kinematic and programmed constraints. The plots are presented using the logarithmic scale. Significantly large violations of the kinematic constraints can be observed for the case of S r L r J nf (Ũ k ¼ 4:39 cm for t ¼ 2:014 s), despite the use of the Baumgarte stabilization method with different parameter settings to eliminate these errors. However, these violations of the programmed constraints are at an acceptable level. One should also pay attention to the different nature of the time courses of violations of the programmed and kinematic constraints. Violations of the programmed constraints, except in the first case, are linear (in the sense of a logarithmic scale), while violations of the kinematic constraints are strongly non-linear. Figure 18 shows the time courses of joint coordinates for the analyzed cases. Examining the calculated time courses, it can be seen that plots obtained for cases with rigid base support are close to each other. Also, the support flexibility impacts significantly courses of the joint coordinates (cases S fl L r J nf , S fl L r J f , S fl L fl J f ). Friction in joints also has noticeable impact on the overall mechanism motion. The flexibility of the coupler has the least influence on motion of individual links of the system. It is interesting because analyzing deformations of the coupler (Fig. 19) it can be concluded that they are large (for t ¼ 2 s deformation of the coupler is around 4 cm for the cases with rigid support). On the other hand, links flexibility has visible impact on forces acting in joints. Figure 20 shows RMS values calculated for the normal and axial forces acting in each joint of the mechanism.
It can be seen that the highest values of the normal and axial forces are obtained for the case in which the coupler (1,3) is flexible. Links flexibility has significant impact on loading of the mechanism and cannot be neglected in numerical simulations.

Conclusions
The general modeling procedure for multibody systems with serial structures of the OLKC and CLKC types subjected to the programmed constraints is presented in the paper. It provides a systematically structured algorithm enabling derivation systems constrained dynamic models with rigid, flexible parts and supports, friction and arbitrary external load. It is based on the joint coordinates formalism and Fig. 17 Constraint violations at the cut-joint R presented in the joint local frame Fig. 18 Time courses of the displacements of the joint coordinates homogeneous transformation matrices. The CoPCoD produces the constrained dynamics, the so called reference dynamics, which enables detailed analysis of the system performance in the desired programmed motion. The underlying dynamics is derived using the GPME algorithm. In simulation tests presented in the paper, for both OLKC and CLKC system models, the influence of the model structure and its properties on the programmed motion were studied. Numerical experiments show that the CoPCoD allows obtaining solutions of the reference dynamics effectively. The computer generation of the model dynamics with the CoPCoD is suitable for various system structures. The presented results indicate how the base flexibility, friction in joints, adoption of rigid or flexible link models impacts motion of the system model. For In our example, the maximum deformations of flexible links are about 10% of the length of links, and their influence on motion of the analyzed systems are negligible. However, the link flexibility has a significant effect on the forces acting in the joints, and cannot be ignored in numerical simulations, even when deformations are small. The future works planned to CoPCoD upgrades to make it more powerful, will be focused on looking for reference dynamics for systems for which the programmed constraints are met with some tolerance, and additional constraints resulting from e.g. design requirements needs to be satisfied.
Funding Open access funding provided by Warsaw University of Technology and University of Bielsko-Biala.
Data availability The C ? ? codes and the datasets generated during the current studies are available from the corresponding author on reasonable request.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.