Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Modelling, simulation and identification of a mobile concrete pump

Modelling, simulation and identification of a mobile concrete pump Mathematical and Computer Modelling of Dynamical Systems, 2015 Vol. 21, No. 2, 180–201, http://dx.doi.org/10.1080/13873954.2014.926277 J. Henikl*, W. Kemmetmüller, M. Bader and A. Kugi Automation and Control Institute, Vienna University of Technology, Vienna, Austria (Received 30 September 2013; accepted 17 May 2014) Due to the light-weight construction of modern large-scale manipulators used, e.g., in mobile concrete pumps, the elasticity of the construction elements plays a significant role in the dynamic behaviour of the system. Therefore, current research is concerned with control strategies for active damping of elastic vibrations and trajectory planning. For this purpose, tailored mathematical models are required. Apart from the mathe- matical modelling, the identification of the model parameters constitutes a challenging task. This is mainly due to the large number of parameters to be identified and, considering the large scale, due to the fact that the boom movement cannot be measured by means of standard sensors. This paper presents a systematic approach for the mathematical modelling and identification of hydraulically actuated large-scale manipulators. The feasibility of the overall approach is demonstrated by means of measurement results of a mobile concrete pump. Keywords: flexible multibody system; hydraulic cylinder; parameter identification; machine vision 1. Introduction New materials and improved construction and production methods enable the design of large-scale manipulators which possess a significantly reduced weight compared to classical designs. Although these new manipulators can carry the same load as traditional constructions, the reduced stiffness makes the systems vulnerable to vibrations. Therefore, the elasticity of the manipulators has to be systematically incorporated into the mathema- tical modelling, the analysis and the controller design. For the modelling of flexible multibody systems, well-developed methods exist in the literature, see e.g., [1,2]. Most investigations in this field are related to flexible robot systems with electromechanical actuators. In large-scale manipulators like mobile concrete pumps, hydraulic actuators comprising hydraulic cylinders and valves are commonly used. Their non- linear dynamic behaviour and the non-linear kinematics of the assembly lead to a complex relationship between piston force and the resulting torque in the joint. Moreover, the dynamics of the pressure in the cylinder chambers is coupled with the piston position and velocity and thus with the motion of the cylinder assembling points on the flexible structure. The combination of flexible multibody systems and hydraulic actuators has been extensively studied in the literature in recent years, see e.g., [3–6]. Therein, different strategies for modelling the elasticity of the boom are proposed. In [4] and [5], the elastic beam elements are approximated by a number of rigid bodies which are connected via *Corresponding author. Email: henikl@acin.tuwien.ac.at This article was originally published with errors. This version has been corrected. Please see Corrigendum http://dx.doi.org/10.1080/13873954.2014.946344 © 2014 Taylor & Francis Mathematical and Computer Modelling of Dynamical Systems 181 joints and spring elements. This leads to a well-known procedure for the derivation of the equations of motion as it is common in robotic applications. The position of the additional joints and the spring rates have, however, no direct relation to real physical parameters. Another approach is to model the elastic segments in form of Euler–Bernoulli beams. This leads to a distributed-parameter representation which has to be approximated by a finite-dimensional model. In this context, the Ritz method is often used in order to obtain a system of ordinary differential equations for the elastic system, see e.g., [7], [3] and [6]. In this work, the latter approach will serve as a basis for the mathematical modelling. In order to ensure fast simulation times of the transient behaviour of the model and to simplify the controller design, the dimension of the model should be kept as small as possible. The right choice of the basis functions for the Ritz method turns out to be the crucial point to achieve a high model accuracy when using only a small number of basis functions. In [7], the eigenfunctions of a homogeneous Euler–Bernoulli beam are chosen as basis functions and [6] recommends the use of Legendre polynomials. These approaches are all feasible, but they do not account for the inhomogeneous structure of the beam elements of the considered application. Therefore, a high number of basis functions are required to accurately describe the beam deflection dynamics. In particular, the static bending of the beam elements cannot be described with sufficient accuracy by a low number of eigenfunctions or polynomial func- tions. In [3], a survey of the choice of the basis functions for hydraulically actuated large-scale manipulators is given. Therein, the use of the static bending line of the homogeneous massless beam with an external unity force is proposed as first basis function. In addition, it is also suggested to utilize the numerically calculated eigenfunctions for inhomogeneous beams. In this paper a set of two basis functions for each beam element is proposed, where the numerically calculated static bending function is used as the first basis function and an orthogonal polynomial serves as a second basis function. As will be shown in the paper, this leads to an adequate approximation of both the static and the dynamic behaviour of the elastic boom. The choice of the elastic degrees of freedom of the boom also has an important influence on the complexity of the resulting equations, which describe the interconnection between the mechanical and the hydraulic part. For example, the model of [4] results in a simple relationship among piston force, joint torque and the rigid body angle while in contrast, the approach of [6] leads to a complex relationship of the generalized forces depending on the elastic deformations. In the present contribution, an adequate formula- tion of the static bending function combined with an orthogonal polynomial ensures a simple kinematic relationship between the motion of the cylinder and the elastic boom. In particular, it is shown that the suggested selection of the degrees of freedom results in a cylinder piston position which solely depends on the rigid body angle. This is of special interest, since this allows a simple incorporation of the static friction of the hydraulic cylinders into the mathematical model. In [8], it is shown that classically designed hydraulic systems used for the actuation of mobile concrete pumps exhibit major weaknesses for the realization of high performance controllers. For this reason, an alternative hydraulic set-up was proposed that fulfils both the dynamic as well as the safety requirements. In this paper, the application of this proposed hydraulic architecture is presumed. This yields a relatively simple mathematical description of the hydraulic components. The identification of the model parameters turns out to be a non-trivial task since the dynamic behaviour of the individual boom segments is highly coupled. A systematic identification procedure is proposed, which is illustrated for the boom of a mobile concrete pump shown in Figure 1. For this purpose, specific static and dynamic experiments were performed, where inclination and pressure sensors, 182 J. Henikl et al. Figure 1. Boom of a mobile concrete pump. strain gauges and incremental encoders were utilized. However, the absolute movement of the boom cannot be measured with this setting. Therefore, the use of non-contact measurement principles like laser, infrared or ultrasonic sensors come into play, but they are practically not feasible since a large number of sensors would be required. This is why we decided to track the movement of particular spots at the boom by means of a camera system by exploiting the methods of machine vision. Commercially available solutions for visual 3D motion like VICON are able to track objects with a precision of 50 μm at 357 Hz with multiple cameras in real time within a closed room as it was used in [9] to control flying quadcopters. Considering, however, the price and the set- up time, such a system would be far oversized, in particular in view of the fact that the motion of the boom for the identification just takes place in a plane of the 3D space. Two high-resolution cameras were used to capture the large working space of the boom. In order to increase the precision, circular markers were attached at the spots to track the boom. This paper is an extended version of [10], where in addition to the derivation of the mathematical model the identification by means of a tailored vision system and the experimental validation is presented in detail. The paper is organized as follows. In Section 2, the derivation of the equations of motion for the elastic boom is shown. The mathematical description of the hydraulic actuation and the friction model are discussed in Section 3.In Section 4, the proposed modelling approach is applied to an industrial mobile concrete pump. For this, the simplifications and assumptions used in the model are summarized, the identification procedure of the model parameters is elaborated and the vision system used to measure the boom movement is described. The simulation results of the mathematical model are compared with measurements in Section 5. Finally, Section 6 gives a short conclusion. 2. Modelling of the mechanical subsystem In this section, the general approach for the derivation of the equations of motion for a large-scale elastic manipulator with multiple beam segments will be presented. Figure 2 illustrates the considered planar flexible manipulator with N ¼ 4 beams and the corre- sponding hydraulic actuators. The degrees of freedom are the rigid-body angles φ and the transversal deflections wðÞ x of the beams, i ¼ 1; ... ; N. i i The overall motion of the system is described with respect to the inertial frame 0 x y . 0 0 0 Furthermore, each beam is equipped with a local coordinate frame 0 x y , where the angle i i i φ describes the rigid body motion of the ith beam with respect to the beam (i−1). By means of a suitable choice of the orientation of the local coordinate frames, the cylinder Mathematical and Computer Modelling of Dynamical Systems 183 x 3 −ϕ 2 −w (x ) 2 2 g x −w (x ) −w (x ) 0 3 3 1 1 3 −ϕ w (x ) 4 4 y 1 0 = 0 0 1 Figure 2. Flexible manipulator with hydraulic actuation. piston position is solely a function of the respective rigid body angle φ . The inputs of the mechanical part are the cylinder forces F . For the derivation of the mathematical model, the kinetic and potential energy of the system are required. For this, the positions and the velocities of each mass particle have to be described in the inertial frame 0 x y . The kinematic relations between the local and 0 0 0 the inertial coordinate frames are given by rotation matrices and translation vectors. The inertial coordinates of a mass particle of the ith beam located at rðÞ x ¼½ x ; wðÞ x , i i i i described in the ith local coordinate frame (subscript), can be calculated by i i i i r ðÞ x ¼ R rðÞ x þ d : (1) i i 0 0 i 0 i 1 2 i The matrix R ¼R R ... R with 0 0 0 i1 cosðÞ φ sinðÞ φ 1 1 1 R ¼ (2) sinðÞ φ cosðÞ φ 1 1 and 1 w L cosðφ Þsinðφ Þ jþ1 j j jþ1 jþ1 R ¼ (3) j 0 w L 1 sinðφ Þ cosðφ Þ j j jþ1 jþ1 for j ¼ 1; ... ; N  1 describes the rotation of the local coordinate frame 0 x y with respect i i i to the inertial frame 0 x y . The lengths of the corresponding beams are denoted by L . 0 0 0 j Similar to [7] small deflections and thus 184 J. Henikl et al. @w @w 0 i i arctan  ¼ wðÞ L (4) @x @x i i x ¼L x ¼L i i i i are presumed. The displacement d between the local coordinate frame 0 x y and the i i i inertial frame 0 x y reads as 0 0 0 i i1 i i1 d ¼ R d þ d ; (5) 0 0 i1 0 T jþ1 with d ¼½ 0; 0 and d ¼ L ; w L for j ¼ 1; ... ; N. j j j Based on the kinematics of the system, the kinetic energy is given by T ¼ T , i¼1 where T is the kinetic energy of the ith beam i i _ _ T ¼ σðÞ x r ðÞ x r ðÞ x dx (6) i i i i i i 0 0 and σðÞ x denotes the position-dependent mass distribution. i i Using the Euler–Bernoulli assumptions, the potential energy due to the beam deflec- tion and gravitation can be described by V ¼ V with i¼1 Z Z L L 2 i i 1 @ wðÞ x i i V ¼ σðÞ x gr ðÞ x dx þ EIðÞ x dx ; (7) i i i i i i i i 0;y 2 @x 0 0 i where g is the acceleration of gravity, E denotes Young’s modulus, IðÞ x is the position- i i i i dependent second moment of inertia and r ðÞ x is the component of r ðÞ x in the i i 0;y 0 direction of y . Viscous damping of the beams is included by means of the Rayleigh dissipation function X i 2 2 R ¼ d φ_ þ d w _ ðÞ x dx ; (8) φi wi i i i i i¼1 with the coefficients d > 0 representing the damping in the joints and d > 0 accounting φ w i i for the damping of the beam deflection movement. Due to the distributed-parameter character of the beams, a mathematical description of the system leads to a set of partial differential equations. To obtain a finite-dimensional approximation, the Ritz method will be employed in the sequel. As already mentioned before, the static bending profiles of the beams are used as first basis functions. The benefit of this choice is the accurate approximation of the elastic deformations at static loads already with the first basis functions alone. Thus, the number of basis functions required for the simulation model can be kept small compared to other choices of the basis functions, e.g., the eigenfunctions of a homogeneous Euler–Bernoulli beam. The shape of the static bending line is determined by the mass distribution of the beam and its external load, which can be represented by a transverse force and torque at the end of the beam. For all beams with the exception of the last one, it can be assumed that the external load is dominating in comparison to the self-weight. Due to the construction of the joints the external load is mainly determined by the transverse forces. Hence, the static bending profiles of the massless beams with position-dependent second moment of inertia IðÞ x i i Mathematical and Computer Modelling of Dynamical Systems 185 and an external force at the free end f serve as basis functions. For the last beam, the bending profile resulting only from its own mass distribution is used. Thus, the static bending profiles w ðÞ x , i ¼ 1; ... ; N  1 are determined by solving i i the differential equations 2 2 @ @ w  ðÞ x i i EIðÞ x ¼ 0 (9) i i 2 2 @x @x i i for an Euler–Bernoulli beam in a left-side clamped configuration, see Figure 3. Due to the assembling of the first hydraulic cylinder at x ¼ x (cf. Figure 2), a supported config- 1 m uration is considered for the first beam, see Figure 4. The boundary conditions for the ordinary differential equations are given by w  ðÞ 0 ¼ 0; (10a) lim w  ðÞ x ¼ 0; (10b) 1 1 x !x 1 m lim w  ðÞ x ¼ 0; (10c) 1 1 x !x þ 1 m @w  ðÞ x @w  ðÞ x 1 1 1 1 lim ¼ lim ; (10d) x !x  x !x þ 1 m @x 1 m @x 1 1 x i w ¯ (x ) i i Figure 3. Static bending profile: beam i ¼ 2; .. . ; N  1. x L m 1 w ¯ (x ) 1 1 Figure 4. Static bending profile: beam 1. 186 J. Henikl et al. @ w  ðÞ x 1 1 EI ðx Þ ¼ 0; (10e) 1 1 x ¼0 @x 2 2 @ w  ðÞ x @ w  ðÞ x 1 1 1 1 lim EI ðÞ x ¼ lim EI ðÞ x ; (10f) 1 1 1 1 2 2 x !x  x !x þ 1 m @x 1 m @x 1 1 @ w  ðÞ x 1 1 EI ðx Þ ¼ 0; (10g) 1 1 x ¼L 1 1 @x @ @ w ðÞ x 1 1 EI ðx Þ ¼ f (10h) 1 1 x ¼L 1 1 L @x @x for the first boom segment and w  ðÞ 0 ¼ 0; (11a) @w  ðÞ x i i ¼ 0; (11b) x ¼0 @x @ w  ðÞ x i i EIðÞ x ¼ 0; (11c) i i x ¼L i i @x @ @ w  ðÞ x i i EIðÞ x ¼ f (11d) i i x ¼L i i L @x @x for the boom segments i= 2, ..., N − 1. The geometric boundary conditions of the joints and the clamped end are represented by Equations (10a–d) and (11a), (11b), respectively. Equations (10e) and (10f) result from the torque-free joints and Equations (10g), (10h) and (11c), (11d) are the boundary conditions for the free ends. For the external force, an arbitrary value f 0 can be presumed. The value of f only scales the result and has no L L influence on the shape of the bending line. The differential equation for the last beam is given by 2 2 @ @ w  ðÞ x N N EI ðÞ x ¼gσ ðÞ x (12) N N N N 2 2 @x @x N N and the boundary conditions, see Figure 5, read as w  ðÞ 0 ¼ 0; (13a) @w  ðÞ x N N ¼ 0; (13b) x ¼0 @x @ w  ðÞ x N N EI ðÞ x ¼ 0; (13c) N N x ¼L N N @x N Mathematical and Computer Modelling of Dynamical Systems 187 w ¯ (x ) N N Figure 5. Static bending profile: beam N. @ @ w  ðÞ x N N EI ðÞ x ¼ 0: (13d) N N x ¼L N N @x @x Since the differential Equations (9)–(13a–d) can only be solved numerically, a polynomial approximation of the static bending profile of each beam is performed, i.e., b b j1 η ðÞ x ¼ ηðÞ x  x x ; (14a) 1 1 m 1 1;j 1 j¼2 b b η ðÞ x ¼ η x ; i ¼ 2; ... ; N ; (14b) i i;j i j¼2 b b b b with the order of the polynomials N and N , respectively. The coefficients η and η 1 i 1;j i;j are determined by minimizing the quadratic error between the numeric solutions of Equations (9)–(13a–d) and the polynomial approximation (14a–b). The formulation of b b η ðÞ x and η ðÞ x has been chosen in such a way that the relative motion of the mounting 1 i 1 i points of the hydraulic actuators on the beams can be solely described by the rigid body angles φ . As a matter of fact, the polynomials (14) fulfil the geometric boundary conditions (10a–d), (11a), (11b) and (13a), (13b). In order to increase the model accuracy especially for the dynamical system behaviour and to account for higher bending modes of the beams, polynomials of the form o 2 o η ðÞ x ¼ðÞ x  x x þ η ðÞ x  x x (15a) 1 1 m 1 m 1 1 1 1;2 o 3 o 2 η ðÞ x ¼ x þ η x ; i ¼ 2; ... ; N (15b) i i i;2 i are used as second basis functions. It can be easily seen that the polynomials (15) also satisfy the geometric boundary conditions. To ensure the orthogonality to the first basis o o functions (14), the coefficients η and η result from the orthogonality condition 1;2 i;2 b o η ðÞ x η ðÞ x dx ¼ 0; i ¼ 1; ... ; N : (16) i i i i i The approximated elastic beam deflections are defined as a superposition of the introduced basis functions b b o o w ¼ h ðÞ t η ðÞ x þ h ðÞ t η ðÞ x ; i ¼ 1; ... ; N ; (17) i i i i i i i 188 J. Henikl et al. b o with the time-varying coefficients h ðÞ t and h ðÞ t . The overall degrees of freedom of i i the mathematical model are composed of the rigid body angles φ and the elastic degrees b o of freedom h , h , which are summarized in the vector i i b b o o q ¼ φ ; ... ; φ ; h ; .. . ; h ; h ; ... ; h : (18) 1 N 1 N 1 N The connection to the hydraulic part is given by the cylinder forces F , which serve as inputs to the mechanical system. Due to the special choice of the local coordinate frames 0 x y , the cylinder piston positions s are solely depending on the rigid body angles φ , i i i p;i i s ¼ fðÞ φ ; i ¼ 1; ... ; N : (19) p;i i The vector of the generalized forces Q can be calculated by means of d’Alembert’s principle in the form @f ðÞ φ @f ðÞ φ 1 1 N N Q ¼ F ; .. . ; F ; 0; ... ; 0; 0; ... ; 0 : (20) 1 N @φ @φ 1 N The functions fðÞ φ are determined by the geometry of the joint construction and are strictly increasing in φ in the operating range of the joints. The application of the Euler–Lagrange equations d @T @T @V @R þ þ ¼ Q (21) _ _ dt @q @q @q @q finally yields the mathematical model of the mechanical part of the system. The resulting equations of motion can be written in the compact form MqðÞq € þcqðÞ ; q_ þ Dq_ þgqðÞ¼ Q; (22) with the positive definite inertia matrix M(q), the vector of Coriolis and centrifugal forces cqðÞ ; q , the damping matrix D and the vector of forces g(q) related to the potential energy. 3. Modelling of the hydraulic subsystem As described in Section 1, large-scale manipulators are typically controlled by means of hydraulic actuators. With the hydraulic architecture proposed in [8], the dynamic behaviour of the hydraulic actuators only depends on the hydraulic cylinders and the proportional directional control valves, see Figure 6. The related differential equations for the chamber pressures p and p take the form 1 2 p_ ¼ A v þ q ; (23a) 1 1 p 1 V þ A s 01 1 p p ¼ A v  q ; (23b) 2 2 p 2 V  A s 02 2 p Mathematical and Computer Modelling of Dynamical Systems 189 A A 1 2 p p 1 2 q q 1 2 p p s t Figure 6. Schematic diagram of the hydraulic actuator. with the bulk modulus of oil β, the piston areas A and A , the offset volumes V and V , 1 2 01 02 and the velocity v ¼ s_ . The volume flows q and q are determined by the valve position p p 1 2 s and the pressure difference, 8 qffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi α A ðÞ s p  p s  0 v1 v s 1 v qffiffi q ¼ (24a) pffiffiffiffiffiffiffiffiffiffiffiffiffiffi α A ðÞ s p  p s < 0; v1 v 1 t v 8 qffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi α A ðÞ s p  p s  0 v2 v 2 t v qffiffi q ¼ (24b) pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi α A ðÞ s p  p s < 0; v2 v s 2 v where α denotes the contraction coefficient, ρ the oil density, p the supply pressure, p the s t tank pressure and A ðÞ s and A ðÞ s are the opening areas of the valve to chamber 1 and v1 v v2 v 2, respectively. The resulting cylinder force F acting on the mechanical part of the system reads as F ¼ p A  p A  F v ; (25) 1 1 2 2 R p where the first part accounts for the pressure force generated by the cylinder and F (v ) R p summarizes the friction forces of the cylinder. It is well known that friction plays an important role in hydraulic actuators. In particular a certain amount of static friction is inevitable for the type of hydraulic cylinders typically used in the considered application. In order to analyse the influence of friction on the dynamic behaviour and the control performance, a static friction model is incorporated into the mathematical model. As can be seen from Figure 7, the friction force is composed of a viscous part r v , a Coulomb V p part r sign(v ) and the static friction r . C p S The static friction phenomenon leads to a switched system behaviour, i.e., if the piston velocity is zero (v = 0) and the absolute value of the sum of the forces on the piston is smaller than r , then the piston sticks and the whole system loses a degree of freedom. S 190 J. Henikl et al. F (v ) R p r v V p −r −r Figure 7. Static friction model. Due to the special choice of the coordinate frame and the basis functions, the detection of zero velocity v ¼ 0 of the ith piston is equivalent to check for φ_ ¼ 0. Furthermore, if p;i i the ith cylinder is sticking, the equations of motion result from Equation (22) by using the constraints φ_ ¼ 0 and φ € ¼ 0. To verify the sticking condition and calculate the absolute i i value of the sum of the forces acting on the piston, the force F of the mechanical part is calculated from the ith row of Equations (20) and (22) N N 1 X X @fðÞ φ € _ _ F ¼ MqðÞ½ i; j q½ j þcqðÞ ; q ½ i þ D½ i; j q½ j þgqðÞ½ i : (26) @φ j¼1 j¼1 Since the functions fðÞ φ are strictly increasing in the operating range, the invertibility of @fðÞ φ =@φ is ensured. Summarizing, the ith piston is sticking if v ¼ 0 and the i p;i i i condition p A  p A  F  r ; (27) 1;i 1;i 2;i 2;i i S;i holds. 4. Parameter identification for a mobile concrete pump In this section, the industrial mobile concrete pump with four joints according to Figure 1 with an operating range of about 40 m is considered. The mathematical model can be derived along the lines of Sections 2 and 3 provided that the following assumptions hold: ● The concrete pump is considered as a planar manipulator. The rotational movement and the lateral beam deflections are neglected. ● The truck’s movement due to the dynamic load of the boom is neglected. ● The boom segments are modelled as Euler–Bernoulli beams. ● The pumping of wet concrete is not considered. However, in view of the last assumption it is worth noting that for the control design the vibrations induced by the pumping are considered as disturbances acting on the system. The model parameters are determined by means of data sheets, design drawings and experimental measurements. In contrast to the parameters of the hydrau- lic components, which are well known from data sheets, the parametrization of the Mathematical and Computer Modelling of Dynamical Systems 191 boom model turns out to be more involved. Besides the geometrical data, the position- dependent mass distribution and the second moment of inertia can be calculated from design drawings. In this context, it is rather difficult to systematically account for externally assembled components on the boom segments like hydraulic cylinders, hydraulic lines and concrete delivery pipes. In particular, their influence on the resulting stiffness of the overall system can hardly be determined analytically. Furthermore, the damping coefficients are unknown and have to be identified by dynamical measurements. 4.1. Experimental set-up In the experimental set-up, in addition to the machine vision system, which will be presented in Section 4.2, the following sensors were used: ● Inclination sensors at both ends of each boom segment. ● Two pressure sensors at each hydraulic cylinder. ● Incremental encoders to measure the cylinder piston movements. Strain gauges attached to the end of each boom segment where the hydraulic actuator is located, in order to determine the beam flection. The strain gauges measure the elongation of the surface of the boom segment. Using the Euler–Bernoulli assumption, the relation between the elongation and the beam flection is given by @ wðÞ x i i ε x ¼ k þ b ; (28) i SG;i SG;i SG;i @x x ¼x i SG;i with the gain k the bias b and the strain gauge position x . In the case of a SG;i SG;i SG;i homogeneous beam, the gain coefficient k is given by the distance of the surface to the SG;i symmetry line of the beam. However, in the case of an inhomogeneous beam it is difficult to determine k analytically. Since in practice it is nearly impossible to attach the strain SG;i gauge to a completely unbent boom segment, the strain gauge signal contains a bias component b , which also depends on the temperature of the beam. Therefore, the gain SG;i and bias coefficients have to be identified by means of measurements. 4.2. Vision system The parameters of the mathematical model of the boom, in particular the mass distribu- tion, the stiffness and damping coefficients and the strain gauge coefficients, could be identified based on static measurements at several configurations of the boom and dynamic experiments with a prescribed excitation of the system. However, due to the large scale of the boom there was the desire to also validate the absolute movement of the boom segments, which cannot be captured by the standard sensors installed. Therefore, a tailored vision system was developed which allows to measure the absolute movement of particular spots on the boom. For this purpose, markers placed at the end of each boom segment and the slewing gear, see Figure 8, were recorded by cameras, and their location in the work space was calculated by means of the respective pixel coordinates on the images. 192 J. Henikl et al. M M M M M 1 2 4 0 3 Figure 8. Markers placed on the boom. In order to capture the large work space of the boom, two 5-mega-pixel greyscale image FireWire cameras were used. The 2452 × 2054 pixel resolution and a frame rate of 15 Hz allow for measurements within a working space of 20 m with an average pixel precision of about 0.01 m/pixel. The open source framework ROS [11] for robotics applications was utilized to record synchronized camera images and to visualize the measurements. After comparing various marker types, including [12–14], circular markers without encoded identifiers were selected. The detection of the circles and their projected elliptical shapes was done by using the standard computer vision library OpenCV. With the algorithm presented in [15] and [12], the precision of the marker position estimation can be increased up to sub-pixel level. This enables the detection of boom movements of even 2 mm. In addition to the markers on the boom, further markers located at fixed positions were used for offline estimation of static camera poses. In order to start the tracking algorithm, all measurement points must be initialized in the camera image. 4.2.1. Feature detection In order to speed up the marker detection process, the algorithm has an initial feature location and predicts the next location throughout the process. For this reason, the marker detection has to be applied only to a small window around the predicted feature location. During the initialization phase, pixel locations and the estimated size of the circular feature have to be defined. The detection algorithm then detects multiple ellipse hypoth- eses in each predicted search window. For this purpose, a function which uses the estimated location and the expected radius sorts all detected ellipses within the search window and places the most likely one first. This ellipse is then used as the basis for the next search window. Figure 9 shows multiple search windows with the most likely ellipse in red. 4.2.2. 3D data generation The geometric relation between the marker position in the workspace and their projected image coordinates are given in this section. The transformation of a point p ¼½p ; p ; p  , defined in the workspace coordinate frame 0 x y z , to the W;x W;y W ;z W W W W W Mathematical and Computer Modelling of Dynamical Systems 193 M M M M M M 2 3 4 0 1 2 (a) (b) Figure 9. Search windows with detected edges in green, ellipses in yellow and the most likely marker ellipse in red. (a) Left camera. (b) Right camera. respective coordinates p ¼½p ; p ; p  in the camera coordinate frame 0 x y z , C;x C;y C;z C C C C 31 which is defined by the camera centre, is given by the translation vector d 2 R and 33 the rotation matrix R 2 R , W W p ¼ R p þ d : (29) C C W C In the literature, this transformation is typically described by the so-called external camera matrix, see e.g., [16]. Assuming a pinhole camera model, the projection of the point p onto the 2D image frame 0 u v is given by I I I "# C;x f þ o p ; u u I u C;z p ¼ ¼ ; (30) I C;y p ; f þ o I v v v C;z with the internal camera parameters f ¼ f =s ; f ¼ f =s ; o and o . This projection is u u v v u v typically referred to as intrinsic camera matrix, see e.g., [16]. The focal length f describes the distance of the camera centre to the sensor chip or the image frame, respectively, s and s represent the effective physical dimensions of one pixel on the camera chip and the distance between the image centre and the origin of the image frame is given by o and o . u v In order to compensate for lens distortions and to be able to use a simple pinhole camera model, every frame captured has to be preprocessed first to get an undistorted image before it is further used in any detection process, see e.g., [16]. The parameters for the compensation of the lens distortions as well as the internal camera parameters were estimated by means of an offline calibration procedure. The calibration was done using a checker-board pattern of known size and the ROS camera calibration node. For the estimation of the three linearly independent quantities parame- W W trizing the rotation matrix R and the three coefficients of the translation vector d ,an C C iterative method based on the Levenberg–Marquardt [17] optimization was utilized, which is implemented in the before mentioned OpenCV library. For this purpose, the workspace coordinates of at least three known reference points p and their visual appearance in the image frame coordinates p are required. Thus, with the transformation (29) and the projection (30), six non-linear equations are given to determine the six unknown coeffi- W W cients for R and d . With the use of additional reference points, the error due to the C C measurement uncertainty of the reference coordinates could be minimized. In Figure 10, five reference points p for j ¼ 1; ... ; 5, used for the estimation of the camera pose, are W ;j shown. The coordinates were measured manually by a laser distance sensor. The reference 194 J. Henikl et al. W,3 W,2 e e p p W,4 W,5 W,1 Figure 10. Reference points. e e points p for l ¼f1; 2; 3; 4g and p for r ¼f1; 2; 3; 5g were used for the left and the W ;l W ;r right camera, respectively. With these results, the relation of the coordinates of a marker centre p described in the workspace coordinate frame and the respective coordinates p is given. The goal here is, however, the inverse problem. A 3D point linked to a 2D image point has to be found. This can be done by creating a straight line, see e.g., [18] that intersects the camera centre and the image frame at the respective 2D image point. Considering Equation (29) and p ¼ 0, the camera centre described in the workspace frame is given by c W W p ¼ R d : (31) W C C With Equation (30) and p ¼ f , the 2D image point at the image frame in front of the C;z m m m camera, defined by the pixel coordinates of the marker centre p ¼½p ; p  ,is I I ;u I ;v located at 2 3 1 m p  o f I;u 6 7 p ¼ 1 m (32) 4 5 p  o f f I;v and m W 1 m W 1 W ^ ¼ðR Þ p ^ ðR Þ d ; (33) W C C C C respectively. With this and Equation (31), the direction of the line is given by 2 3 1 m p  o I;u 6 7 m c W 6 7 v ¼ p  p ¼ f R 1 m ¼ f v (34) W W W W C 4 5 p  o f I;v and thus, the line could be described by the equation Mathematical and Computer Modelling of Dynamical Systems 195 l c p ðλÞ¼ p þ λv (35) W W for λ 2 R. The final 3D marker location p can be computed under the assumption that the boom and the marker moves only in the plane defined by n p ¼ λ ; (36) with the normal vector n to the plane and a scalar λ defining the position of the plane. The normal vector n and the scalar λ were determined by means of the reference points lying in the plane, see Figure 10. With this, the scalar λ has to be found that defines the c m intersection of the line (35) with the plane (36), p þ λ v ¼ p . Considering τ W W W T c n ðp þ λ v Þ¼ λ ; λ is given by τ W  τ T c λ  n p λ ¼ ; (37) n v and therefore p reads as T c T c λ  n p λ  n p m c  W c  W p ¼ p þ v ¼ p þ v : (38) W W W W W T T n v n v W W 4.3. Identification procedure Due to the model complexity and the high number of system parameters, the identification has to be performed in several steps. The following procedure was applied to the simulation model of the mobile concrete pump. First, a fine tuning of the mass distribution was obtained by means of static experi- ments. For this purpose, the cylinder forces respectively the chamber pressures and the inclination at the beginning of each boom segment were measured at different angles φ ¼ φ  of the first joint and a sprawled boom, i.e., φ ¼ φ ¼ φ ¼ 0. For the pre- 1;m 1 2 3 4 scribed constant angles φ  , the corresponding measured cylinder forces and inclinations 1;m are referred to as F ðφ  Þ and ψ  φ  ; i ¼ 1; ... ; N. In order to compare the i;m 1;m b;i;m 1;m measurement results with the simulation model, the stationary solution of Equations (20) and (22) is determined, i.e., gðq Þ¼ Q and, given φ ¼ φ ¼ φ ¼ 0, solved for 2 3 4 b b o o h ; ... ; h ; h ; ... ; h and F ; ... ; F as a function of φ  . In general, the inclination at the 1 N 1 1 N 1 N beginning of the boom segment reads as i i1 X X ψ ¼ φ þ w ðL Þ: (39) n n b;i n n¼1 n¼1 Thus, for φ ¼ φ ¼ φ ¼ 0 and w (x ) according to Equation (17), the stationary 2 3 4 n n relation of Equation (39) takes the form i1 b b o o ψ ¼ φ þ h η ðÞ L þ h η ðÞ L : (40) b;i 1 n n n n n n n¼1 196 J. Henikl et al. With Equation (40), the stationary cylinder forces F ; i ¼ 1; ... ; N, can be parame- trized in terms of the respective inclination, i.e., F ðψ  Þ. Now the nominal mass distribu- i b;i tion of the boom segments σ ðx Þ; i ¼ 1; ... ; N is scaled with a multiplicative constant i i ; i ¼ 1; ... ; N, which serves as a parameter to be identified. For this, the least squares problem min F ðψ  Þ F ðψ  Þ (41) i i;m b;i;m b;i;m has to be recursively solved beginning at the last boom segment i = N. In order to identify the stiffness and the damping behaviour of the boom, dynamic experiments with a defined excitation of the system were performed in the next step. With the measurements of the cylinder forces, the strain gauges and the marker movements of the machine vision system, the mean equivalent Young’s modulus E as well as the damping coefficients d and d were determined such that the simulation model matches φ w the occurring oscillation frequencies and the decay of the envelopes. The strain gauge coefficients k and b were also identified by means of static measurements, where a SG;i SG;i similar procedure as for the adjustment of the mass distribution was employed. 5. Experimental results The following subsection presents the validation of the simulation model by means of experimental results. The numerically calculated basis functions (14) and (15), employed in the simulation model, are illustrated in Figure 11.In Figure 12, the results of static measurements in comparison to the simulation model are given. The cylinder forces and beam flections have been measured at different angles φ of the first joint and a sprawled boom with joint angles φ , φ and φ close to 0. The cylinder forces are normalized to the 2 3 4 1 1 0 0 x = x 1 m −1 −1 0 0.2 0.4 0.6 0.8 1 00.2 0.4 0.6 0.8 1 x /L x /L 1 1 2 2 1 1 0 0 −1 −1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x /L x /L 3 3 4 4 1st basis function 2nd basis function Figure 11. Basis functions of the boom segments normalized to the beam length and to the maximum function value. η (x )/max|η (x )| η (x )/max|η (x )| 3 3 3 3 1 1 1 1 η (x )/max|η (x )| η (x )/max|η (x )| 4 4 4 4 2 2 2 2 Mathematical and Computer Modelling of Dynamical Systems 197 0. 6 0.5 0. 4 0.3 0. 2 0.1 0. 5 0.3 0. 3 0.2 0.1 0. 1 500 400 250 150 020 40 60 80 020 40 60 80 ϕ [°] ϕ [°] 1 1 Simulation Measurement Figure 12. Static measurements at a sprawled boom: cylinder forces normalized to the maximum forces and strain gauge measurements as a function of the first joint angle. maximum forces, which are determined by the supply pressure and the piston area, F ¼ p A . i;N s 1:i In order to validate the dynamic system behaviour, the system was excited by a sudden release of a 90 kg load at the end of the boom. In Figure 13, the comparison of the measured and simulated response of the normalized cylinder forces, the strain gauge signals and relative movements of the cylinder pistons at a typical boom configuration on the construction site is shown, see Figure 1. The stationary deviations of the strain gauge signals can be explained by thermal expansion due to different experiment days. Finally, in Figure 14 the measurement results of the machine vision system are presented. The simulated and tracked movement of the markers, described in the inertial frame, is illustrated. It can be seen that the amplitudes of the oscillations agree very well. With this, an essential requirement for the mathematical model, which serves as a basis for the design of active damping control strategies, is fulfilled. Looking at the absolute deviations to the inertial frame, one has to bear in mind that already small deviations of the angle w ðL Þ at the end of the first boom segment have a big influence on the absolute marker positions. Furthermore, due to lens distortions at the outer image section as well as the restricted validity of the assumption that all markers move in the same plane, the precision of the absolute machine vision measurements is limited. However, considering the size of the boom, the biggest vertical and horizontal deviations of about 60 cm for y , m 4 and 15 cm for x can be considered small in relation to the horizontal range of 27 m. m;3 ε [μm/m] ε [μm/m] F /F F /F 3 3 3,N 1 1,N ε [μm/m] ε [μm/m] F /F F /F 4 2 4 4,N 2 2,N 198 J. Henikl et al. 0.62 0.56 0.58 0.52 0.48 0.54 0.4 0.44 0.4 0.3 0.36 0.2 –50 0.6 0.3 0.4 0.2 0.1 0.2 0.6 0.8 0.6 0.4 0.4 0.2 0.2 05 10 15 20 0 510 15 20 t [s] t [s] Simulation Measurement Figure 13. Load release at the end-effector: cylinder forces, strain gauge signals and cylinder piston movements. Since in view of the envisaged controller design the main goal of the model is to capture the essential dynamic behaviour of the mobile concrete pump, the impreciseness in terms of the absolute deviations is of minor importance. Figure 14 shows the relatively high resolution and low noise level of the signals, taking into account the large work space. Figure 15 depicts the relative vertical movement of a marker attached to the slewing gear. It can be seen that the vision system is even able to capture the 2 mm amplitude of the reaction of the truck caused by the sudden release of the load. With the use of the effective ellipse search algorithm, the image processing time for each frame is far lower than the sampling time defined by the frame rate. Hence, it would even be possible to use this system for real-time applications. Δs [mm] Δs [mm] ε [μm/m] ε [μm/m] F /F F /F p,3 p,1 3 1 3 3,N 1 1,N Δs [mm] Δs [mm] ε [μm/m] ε [μm/m] F /F F /F p,4 p,2 4 2 4 4,N 2 2,N Mathematical and Computer Modelling of Dynamical Systems 199 14.7 6.77 6.75 14.5 6.73 14.3 27.5 21.6 21.55 27.3 21.5 21.45 27.1 4.1 5.3 5.2 4.05 5.1 2.6 −2.4 2.4 −2.8 2.2 −3.2 05 10 15 20 05 10 15 20 t [s] t [s] Simulation Left camera Right camera Figure 14. Load release at the end-effector: horizontal and vertical distance of the markers to the inertial frame. −1 −2 −3 02468 10 12 14 16 18 20 t [s] Figure 15. Load release at the end-effector: relative vertical movement of a marker attached to the truck. 6. Conclusions In this paper, a systematic approach for the mathematical modelling of large-scale light- weight manipulators with hydraulic actuation was proposed. It was shown that a suitable choice of the basis functions for the approximation of the elastic bending of the beams has significant advantages for the incorporation of static friction into the model. Furthermore, simulation and measurement results of a mobile concrete pump y [m] y [m] x [m] x [m] m,3 m,1 m,3 m,1 Δ [mm] ym,0 y [m] y [m] x [m] x [m] m,4 m,2 m,4 m,2 200 J. Henikl et al. show the accuracy and the feasibility of the proposed modelling approach already for a small number of basis functions. A procedure for the identification of the uncertain model parameters is given. For this purpose, a tailored machine vision system was presented for the measurement of the absolute boom movement. It was demonstrated that this system enables the tracking of motions with amplitudes of approximately 2 mm even for a large observation area. Acknowledgements The authors want to express their gratitude to Helmut Brunner for the implementation of the vision system algorithm as well as the recording and evaluation of the camera images. Funding Part of the research of this work was supported by the company TTControl and the Austrian Research Promotion Agency (FFG) within the project MOBILE [Project no. 2426433]. Supplemental data Supplemental data for this article can be accessed here. Notes 1. VICON: http://www.vicon.com 2. For the sake of readability, the time dependence of the degrees of freedom is not explicitly stated in this paper. 3. Allide Vision: Pike F-505, http://www.alliedvisiontec.com 4. ROS (Robot Operating System): http://www.ros.org 5. OpenCV (Open Source Computer Vision): http://opencv.willowgarage.com/wiki/ 6. ROS Camera Calibration: www.ros.org/wiki/camera_calibration 6 1 7. With the linear thermal expansion coefficient for steel of α  11 ... 13  10 K a T ;Steel temperature difference over night of 10 K yields an expansion of about 110 ... 130 μm=m: 8. A video illustrating the capturing of the vertical motion of the markers is available on http:// www.acin.tuwien.ac.at/fileadmin/cds/videos/Motion_Capturing_Mobile_Concrete_Pump.mp4 References [1] A.A. Shabana, Dynamics of Multibody Systems, Cambridge University Press, New York, 2005. [2] H. Bremer, Elastic Multibody Dynamics: A Direct Ritz Approach, Springer, Heidelberg, 2008. [3] M. Schneider, Modellbildung, Simulation und nichtlineare Regelung elastischer, hydraulisch angetriebener Großraummanipulatoren, VDI Fortschrittsberichte Reihe 8, VDI Verlag, Düsseldorf, 1999. [4] H. Aschemann, O. Sawodny, A. Bulach, and H.P. Hofer, Model based trajectory control of a flexible turntable ladder, Proceedings of the American Control Conference Vol. 2, Anchorage, AK, 8–10 May 2002, pp. 921–926. [5] S. Lambeck, O. Sawodny, and E. Arnold, Trajectory tracking control for a new generation of fire rescue turntable ladders, Proceedings of the 2nd IEEE Conference on Robotics, Automation and Mechatronics, Bangkok, 1–3 June 2006, pp. 847–852. [6] M. Dibold, J. Gerstmayr, R. Stadlmayr, H. Irschik, and K. Schlacher, Dynamics of multibody systems including hydraulic actuators and feedback control, Proceedings of the 3rd European Conference on Structural Control Vol. 2, Vienna, 12–15 July 2004, pp. 39–42. [7] A. De Luca and B. Siciliano, Closed-form dynamic model of planar multilink lightweight robots, IEEE T Syst. Man Cy. 21 (1991), pp. 826–839. doi:10.1109/21.108300. [8] J. Henikl, W. Kemmetmüller, and A. Kugi, Modeling and control of a mobile concrete pump, Proceedings of the 6th IFAC Symposium of Mechatronic Systems, Hangzhou, 10–12 April 2013, pp. 91–98. Mathematical and Computer Modelling of Dynamical Systems 201 [9] R. Mahony, V. Kumar, and P. Corke, Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor, IEEE Robot. Autom. Mag. 19 (2012), pp. 20–32. doi:10.1109/ MRA.2012.2206474. [10] J. Henikl, W. Kemmetmuüller, and A. Kugi, Modeling and simulation of large-scale manip- ulators with hydraulic actuation, Proceedings of the 7th Vienna Conference on Mathematical Modeling, Vienna, 14–17 February 2012, pp. 780–785. [11] M. Quigley, K. Conley, B. Gerkey, J. Faust, T. Foote, J. Leibs, R. Wheeler, and A.Y. Ng, ROS: An open-source robot operating system, ICRA Workshop on Open Source Software, Kobe, 12–17 May 2009. [12] X. Zhang, S. Fronz, and N. Navab, Visual marker detection and decoding in AR systems: A comparative study, Proceedings of the 1st International Symposium on Mixed and Augmented Reality, Darmstadt, 30 September–1 October 2002, pp. 97–106. [13] L. Naimark and E. Foxlin, Circular data matrix fiducial system and robust image processing for a wearable vision-inertial self-tracker, Proceedings of the 1st International Symposium on Mixed and Augmented Reality, Darmstadt, 30 September–1 October 2002, pp. 27–36. [14] F. Ababsa and M. Mallem, A robust circular fiducial detection technique and real-time 3D camera tracking, J. Multimedia. 3 (2008), pp. 34–41. doi:10.4304/jmm.3.4.34-41. [15] M. Zillich and J. Matas, Ellipse detection using efficient grouping of arc segments, Proceedings of the 27th Workshop of the Austrian Association for Pattern Recognition (OAGM/AAPR), Laxenburg, 5–6 June 2003, pp. 143–148. [16] E. Trucco and A. Verri, Introductory Techniques for 3-D Computer Vision, Prentice Hall PTR, Upper Saddle River, NJ, 1998. [17] D. Grest, T. Petersen, and V. Krüger, A comparison of iterative 2D-3D pose estimation methods for real-time applications,in Image Analysis, A.B. Salberg, J. Hardeberg, and R. Jenssen, eds., Springer, Heidelberg, 2009, pp. 706–715. [18] C. Wohler, 3D Computer Vision: Efficient Methods and Applications, Springer, Heidelberg, http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematical and Computer Modelling of Dynamical Systems Taylor & Francis

Modelling, simulation and identification of a mobile concrete pump

Download PDF
22 pages

Loading next page...
 
/lp/taylor-francis/modelling-simulation-and-identification-of-a-mobile-concrete-pump-c1gUJF7gZM

References (20)

Publisher
Taylor & Francis
Copyright
© 2014 Taylor & Francis
ISSN
1744-5051
eISSN
1387-3954
DOI
10.1080/13873954.2014.926277
Publisher site
See Article on Publisher Site

Abstract

Mathematical and Computer Modelling of Dynamical Systems, 2015 Vol. 21, No. 2, 180–201, http://dx.doi.org/10.1080/13873954.2014.926277 J. Henikl*, W. Kemmetmüller, M. Bader and A. Kugi Automation and Control Institute, Vienna University of Technology, Vienna, Austria (Received 30 September 2013; accepted 17 May 2014) Due to the light-weight construction of modern large-scale manipulators used, e.g., in mobile concrete pumps, the elasticity of the construction elements plays a significant role in the dynamic behaviour of the system. Therefore, current research is concerned with control strategies for active damping of elastic vibrations and trajectory planning. For this purpose, tailored mathematical models are required. Apart from the mathe- matical modelling, the identification of the model parameters constitutes a challenging task. This is mainly due to the large number of parameters to be identified and, considering the large scale, due to the fact that the boom movement cannot be measured by means of standard sensors. This paper presents a systematic approach for the mathematical modelling and identification of hydraulically actuated large-scale manipulators. The feasibility of the overall approach is demonstrated by means of measurement results of a mobile concrete pump. Keywords: flexible multibody system; hydraulic cylinder; parameter identification; machine vision 1. Introduction New materials and improved construction and production methods enable the design of large-scale manipulators which possess a significantly reduced weight compared to classical designs. Although these new manipulators can carry the same load as traditional constructions, the reduced stiffness makes the systems vulnerable to vibrations. Therefore, the elasticity of the manipulators has to be systematically incorporated into the mathema- tical modelling, the analysis and the controller design. For the modelling of flexible multibody systems, well-developed methods exist in the literature, see e.g., [1,2]. Most investigations in this field are related to flexible robot systems with electromechanical actuators. In large-scale manipulators like mobile concrete pumps, hydraulic actuators comprising hydraulic cylinders and valves are commonly used. Their non- linear dynamic behaviour and the non-linear kinematics of the assembly lead to a complex relationship between piston force and the resulting torque in the joint. Moreover, the dynamics of the pressure in the cylinder chambers is coupled with the piston position and velocity and thus with the motion of the cylinder assembling points on the flexible structure. The combination of flexible multibody systems and hydraulic actuators has been extensively studied in the literature in recent years, see e.g., [3–6]. Therein, different strategies for modelling the elasticity of the boom are proposed. In [4] and [5], the elastic beam elements are approximated by a number of rigid bodies which are connected via *Corresponding author. Email: henikl@acin.tuwien.ac.at This article was originally published with errors. This version has been corrected. Please see Corrigendum http://dx.doi.org/10.1080/13873954.2014.946344 © 2014 Taylor & Francis Mathematical and Computer Modelling of Dynamical Systems 181 joints and spring elements. This leads to a well-known procedure for the derivation of the equations of motion as it is common in robotic applications. The position of the additional joints and the spring rates have, however, no direct relation to real physical parameters. Another approach is to model the elastic segments in form of Euler–Bernoulli beams. This leads to a distributed-parameter representation which has to be approximated by a finite-dimensional model. In this context, the Ritz method is often used in order to obtain a system of ordinary differential equations for the elastic system, see e.g., [7], [3] and [6]. In this work, the latter approach will serve as a basis for the mathematical modelling. In order to ensure fast simulation times of the transient behaviour of the model and to simplify the controller design, the dimension of the model should be kept as small as possible. The right choice of the basis functions for the Ritz method turns out to be the crucial point to achieve a high model accuracy when using only a small number of basis functions. In [7], the eigenfunctions of a homogeneous Euler–Bernoulli beam are chosen as basis functions and [6] recommends the use of Legendre polynomials. These approaches are all feasible, but they do not account for the inhomogeneous structure of the beam elements of the considered application. Therefore, a high number of basis functions are required to accurately describe the beam deflection dynamics. In particular, the static bending of the beam elements cannot be described with sufficient accuracy by a low number of eigenfunctions or polynomial func- tions. In [3], a survey of the choice of the basis functions for hydraulically actuated large-scale manipulators is given. Therein, the use of the static bending line of the homogeneous massless beam with an external unity force is proposed as first basis function. In addition, it is also suggested to utilize the numerically calculated eigenfunctions for inhomogeneous beams. In this paper a set of two basis functions for each beam element is proposed, where the numerically calculated static bending function is used as the first basis function and an orthogonal polynomial serves as a second basis function. As will be shown in the paper, this leads to an adequate approximation of both the static and the dynamic behaviour of the elastic boom. The choice of the elastic degrees of freedom of the boom also has an important influence on the complexity of the resulting equations, which describe the interconnection between the mechanical and the hydraulic part. For example, the model of [4] results in a simple relationship among piston force, joint torque and the rigid body angle while in contrast, the approach of [6] leads to a complex relationship of the generalized forces depending on the elastic deformations. In the present contribution, an adequate formula- tion of the static bending function combined with an orthogonal polynomial ensures a simple kinematic relationship between the motion of the cylinder and the elastic boom. In particular, it is shown that the suggested selection of the degrees of freedom results in a cylinder piston position which solely depends on the rigid body angle. This is of special interest, since this allows a simple incorporation of the static friction of the hydraulic cylinders into the mathematical model. In [8], it is shown that classically designed hydraulic systems used for the actuation of mobile concrete pumps exhibit major weaknesses for the realization of high performance controllers. For this reason, an alternative hydraulic set-up was proposed that fulfils both the dynamic as well as the safety requirements. In this paper, the application of this proposed hydraulic architecture is presumed. This yields a relatively simple mathematical description of the hydraulic components. The identification of the model parameters turns out to be a non-trivial task since the dynamic behaviour of the individual boom segments is highly coupled. A systematic identification procedure is proposed, which is illustrated for the boom of a mobile concrete pump shown in Figure 1. For this purpose, specific static and dynamic experiments were performed, where inclination and pressure sensors, 182 J. Henikl et al. Figure 1. Boom of a mobile concrete pump. strain gauges and incremental encoders were utilized. However, the absolute movement of the boom cannot be measured with this setting. Therefore, the use of non-contact measurement principles like laser, infrared or ultrasonic sensors come into play, but they are practically not feasible since a large number of sensors would be required. This is why we decided to track the movement of particular spots at the boom by means of a camera system by exploiting the methods of machine vision. Commercially available solutions for visual 3D motion like VICON are able to track objects with a precision of 50 μm at 357 Hz with multiple cameras in real time within a closed room as it was used in [9] to control flying quadcopters. Considering, however, the price and the set- up time, such a system would be far oversized, in particular in view of the fact that the motion of the boom for the identification just takes place in a plane of the 3D space. Two high-resolution cameras were used to capture the large working space of the boom. In order to increase the precision, circular markers were attached at the spots to track the boom. This paper is an extended version of [10], where in addition to the derivation of the mathematical model the identification by means of a tailored vision system and the experimental validation is presented in detail. The paper is organized as follows. In Section 2, the derivation of the equations of motion for the elastic boom is shown. The mathematical description of the hydraulic actuation and the friction model are discussed in Section 3.In Section 4, the proposed modelling approach is applied to an industrial mobile concrete pump. For this, the simplifications and assumptions used in the model are summarized, the identification procedure of the model parameters is elaborated and the vision system used to measure the boom movement is described. The simulation results of the mathematical model are compared with measurements in Section 5. Finally, Section 6 gives a short conclusion. 2. Modelling of the mechanical subsystem In this section, the general approach for the derivation of the equations of motion for a large-scale elastic manipulator with multiple beam segments will be presented. Figure 2 illustrates the considered planar flexible manipulator with N ¼ 4 beams and the corre- sponding hydraulic actuators. The degrees of freedom are the rigid-body angles φ and the transversal deflections wðÞ x of the beams, i ¼ 1; ... ; N. i i The overall motion of the system is described with respect to the inertial frame 0 x y . 0 0 0 Furthermore, each beam is equipped with a local coordinate frame 0 x y , where the angle i i i φ describes the rigid body motion of the ith beam with respect to the beam (i−1). By means of a suitable choice of the orientation of the local coordinate frames, the cylinder Mathematical and Computer Modelling of Dynamical Systems 183 x 3 −ϕ 2 −w (x ) 2 2 g x −w (x ) −w (x ) 0 3 3 1 1 3 −ϕ w (x ) 4 4 y 1 0 = 0 0 1 Figure 2. Flexible manipulator with hydraulic actuation. piston position is solely a function of the respective rigid body angle φ . The inputs of the mechanical part are the cylinder forces F . For the derivation of the mathematical model, the kinetic and potential energy of the system are required. For this, the positions and the velocities of each mass particle have to be described in the inertial frame 0 x y . The kinematic relations between the local and 0 0 0 the inertial coordinate frames are given by rotation matrices and translation vectors. The inertial coordinates of a mass particle of the ith beam located at rðÞ x ¼½ x ; wðÞ x , i i i i described in the ith local coordinate frame (subscript), can be calculated by i i i i r ðÞ x ¼ R rðÞ x þ d : (1) i i 0 0 i 0 i 1 2 i The matrix R ¼R R ... R with 0 0 0 i1 cosðÞ φ sinðÞ φ 1 1 1 R ¼ (2) sinðÞ φ cosðÞ φ 1 1 and 1 w L cosðφ Þsinðφ Þ jþ1 j j jþ1 jþ1 R ¼ (3) j 0 w L 1 sinðφ Þ cosðφ Þ j j jþ1 jþ1 for j ¼ 1; ... ; N  1 describes the rotation of the local coordinate frame 0 x y with respect i i i to the inertial frame 0 x y . The lengths of the corresponding beams are denoted by L . 0 0 0 j Similar to [7] small deflections and thus 184 J. Henikl et al. @w @w 0 i i arctan  ¼ wðÞ L (4) @x @x i i x ¼L x ¼L i i i i are presumed. The displacement d between the local coordinate frame 0 x y and the i i i inertial frame 0 x y reads as 0 0 0 i i1 i i1 d ¼ R d þ d ; (5) 0 0 i1 0 T jþ1 with d ¼½ 0; 0 and d ¼ L ; w L for j ¼ 1; ... ; N. j j j Based on the kinematics of the system, the kinetic energy is given by T ¼ T , i¼1 where T is the kinetic energy of the ith beam i i _ _ T ¼ σðÞ x r ðÞ x r ðÞ x dx (6) i i i i i i 0 0 and σðÞ x denotes the position-dependent mass distribution. i i Using the Euler–Bernoulli assumptions, the potential energy due to the beam deflec- tion and gravitation can be described by V ¼ V with i¼1 Z Z L L 2 i i 1 @ wðÞ x i i V ¼ σðÞ x gr ðÞ x dx þ EIðÞ x dx ; (7) i i i i i i i i 0;y 2 @x 0 0 i where g is the acceleration of gravity, E denotes Young’s modulus, IðÞ x is the position- i i i i dependent second moment of inertia and r ðÞ x is the component of r ðÞ x in the i i 0;y 0 direction of y . Viscous damping of the beams is included by means of the Rayleigh dissipation function X i 2 2 R ¼ d φ_ þ d w _ ðÞ x dx ; (8) φi wi i i i i i¼1 with the coefficients d > 0 representing the damping in the joints and d > 0 accounting φ w i i for the damping of the beam deflection movement. Due to the distributed-parameter character of the beams, a mathematical description of the system leads to a set of partial differential equations. To obtain a finite-dimensional approximation, the Ritz method will be employed in the sequel. As already mentioned before, the static bending profiles of the beams are used as first basis functions. The benefit of this choice is the accurate approximation of the elastic deformations at static loads already with the first basis functions alone. Thus, the number of basis functions required for the simulation model can be kept small compared to other choices of the basis functions, e.g., the eigenfunctions of a homogeneous Euler–Bernoulli beam. The shape of the static bending line is determined by the mass distribution of the beam and its external load, which can be represented by a transverse force and torque at the end of the beam. For all beams with the exception of the last one, it can be assumed that the external load is dominating in comparison to the self-weight. Due to the construction of the joints the external load is mainly determined by the transverse forces. Hence, the static bending profiles of the massless beams with position-dependent second moment of inertia IðÞ x i i Mathematical and Computer Modelling of Dynamical Systems 185 and an external force at the free end f serve as basis functions. For the last beam, the bending profile resulting only from its own mass distribution is used. Thus, the static bending profiles w ðÞ x , i ¼ 1; ... ; N  1 are determined by solving i i the differential equations 2 2 @ @ w  ðÞ x i i EIðÞ x ¼ 0 (9) i i 2 2 @x @x i i for an Euler–Bernoulli beam in a left-side clamped configuration, see Figure 3. Due to the assembling of the first hydraulic cylinder at x ¼ x (cf. Figure 2), a supported config- 1 m uration is considered for the first beam, see Figure 4. The boundary conditions for the ordinary differential equations are given by w  ðÞ 0 ¼ 0; (10a) lim w  ðÞ x ¼ 0; (10b) 1 1 x !x 1 m lim w  ðÞ x ¼ 0; (10c) 1 1 x !x þ 1 m @w  ðÞ x @w  ðÞ x 1 1 1 1 lim ¼ lim ; (10d) x !x  x !x þ 1 m @x 1 m @x 1 1 x i w ¯ (x ) i i Figure 3. Static bending profile: beam i ¼ 2; .. . ; N  1. x L m 1 w ¯ (x ) 1 1 Figure 4. Static bending profile: beam 1. 186 J. Henikl et al. @ w  ðÞ x 1 1 EI ðx Þ ¼ 0; (10e) 1 1 x ¼0 @x 2 2 @ w  ðÞ x @ w  ðÞ x 1 1 1 1 lim EI ðÞ x ¼ lim EI ðÞ x ; (10f) 1 1 1 1 2 2 x !x  x !x þ 1 m @x 1 m @x 1 1 @ w  ðÞ x 1 1 EI ðx Þ ¼ 0; (10g) 1 1 x ¼L 1 1 @x @ @ w ðÞ x 1 1 EI ðx Þ ¼ f (10h) 1 1 x ¼L 1 1 L @x @x for the first boom segment and w  ðÞ 0 ¼ 0; (11a) @w  ðÞ x i i ¼ 0; (11b) x ¼0 @x @ w  ðÞ x i i EIðÞ x ¼ 0; (11c) i i x ¼L i i @x @ @ w  ðÞ x i i EIðÞ x ¼ f (11d) i i x ¼L i i L @x @x for the boom segments i= 2, ..., N − 1. The geometric boundary conditions of the joints and the clamped end are represented by Equations (10a–d) and (11a), (11b), respectively. Equations (10e) and (10f) result from the torque-free joints and Equations (10g), (10h) and (11c), (11d) are the boundary conditions for the free ends. For the external force, an arbitrary value f 0 can be presumed. The value of f only scales the result and has no L L influence on the shape of the bending line. The differential equation for the last beam is given by 2 2 @ @ w  ðÞ x N N EI ðÞ x ¼gσ ðÞ x (12) N N N N 2 2 @x @x N N and the boundary conditions, see Figure 5, read as w  ðÞ 0 ¼ 0; (13a) @w  ðÞ x N N ¼ 0; (13b) x ¼0 @x @ w  ðÞ x N N EI ðÞ x ¼ 0; (13c) N N x ¼L N N @x N Mathematical and Computer Modelling of Dynamical Systems 187 w ¯ (x ) N N Figure 5. Static bending profile: beam N. @ @ w  ðÞ x N N EI ðÞ x ¼ 0: (13d) N N x ¼L N N @x @x Since the differential Equations (9)–(13a–d) can only be solved numerically, a polynomial approximation of the static bending profile of each beam is performed, i.e., b b j1 η ðÞ x ¼ ηðÞ x  x x ; (14a) 1 1 m 1 1;j 1 j¼2 b b η ðÞ x ¼ η x ; i ¼ 2; ... ; N ; (14b) i i;j i j¼2 b b b b with the order of the polynomials N and N , respectively. The coefficients η and η 1 i 1;j i;j are determined by minimizing the quadratic error between the numeric solutions of Equations (9)–(13a–d) and the polynomial approximation (14a–b). The formulation of b b η ðÞ x and η ðÞ x has been chosen in such a way that the relative motion of the mounting 1 i 1 i points of the hydraulic actuators on the beams can be solely described by the rigid body angles φ . As a matter of fact, the polynomials (14) fulfil the geometric boundary conditions (10a–d), (11a), (11b) and (13a), (13b). In order to increase the model accuracy especially for the dynamical system behaviour and to account for higher bending modes of the beams, polynomials of the form o 2 o η ðÞ x ¼ðÞ x  x x þ η ðÞ x  x x (15a) 1 1 m 1 m 1 1 1 1;2 o 3 o 2 η ðÞ x ¼ x þ η x ; i ¼ 2; ... ; N (15b) i i i;2 i are used as second basis functions. It can be easily seen that the polynomials (15) also satisfy the geometric boundary conditions. To ensure the orthogonality to the first basis o o functions (14), the coefficients η and η result from the orthogonality condition 1;2 i;2 b o η ðÞ x η ðÞ x dx ¼ 0; i ¼ 1; ... ; N : (16) i i i i i The approximated elastic beam deflections are defined as a superposition of the introduced basis functions b b o o w ¼ h ðÞ t η ðÞ x þ h ðÞ t η ðÞ x ; i ¼ 1; ... ; N ; (17) i i i i i i i 188 J. Henikl et al. b o with the time-varying coefficients h ðÞ t and h ðÞ t . The overall degrees of freedom of i i the mathematical model are composed of the rigid body angles φ and the elastic degrees b o of freedom h , h , which are summarized in the vector i i b b o o q ¼ φ ; ... ; φ ; h ; .. . ; h ; h ; ... ; h : (18) 1 N 1 N 1 N The connection to the hydraulic part is given by the cylinder forces F , which serve as inputs to the mechanical system. Due to the special choice of the local coordinate frames 0 x y , the cylinder piston positions s are solely depending on the rigid body angles φ , i i i p;i i s ¼ fðÞ φ ; i ¼ 1; ... ; N : (19) p;i i The vector of the generalized forces Q can be calculated by means of d’Alembert’s principle in the form @f ðÞ φ @f ðÞ φ 1 1 N N Q ¼ F ; .. . ; F ; 0; ... ; 0; 0; ... ; 0 : (20) 1 N @φ @φ 1 N The functions fðÞ φ are determined by the geometry of the joint construction and are strictly increasing in φ in the operating range of the joints. The application of the Euler–Lagrange equations d @T @T @V @R þ þ ¼ Q (21) _ _ dt @q @q @q @q finally yields the mathematical model of the mechanical part of the system. The resulting equations of motion can be written in the compact form MqðÞq € þcqðÞ ; q_ þ Dq_ þgqðÞ¼ Q; (22) with the positive definite inertia matrix M(q), the vector of Coriolis and centrifugal forces cqðÞ ; q , the damping matrix D and the vector of forces g(q) related to the potential energy. 3. Modelling of the hydraulic subsystem As described in Section 1, large-scale manipulators are typically controlled by means of hydraulic actuators. With the hydraulic architecture proposed in [8], the dynamic behaviour of the hydraulic actuators only depends on the hydraulic cylinders and the proportional directional control valves, see Figure 6. The related differential equations for the chamber pressures p and p take the form 1 2 p_ ¼ A v þ q ; (23a) 1 1 p 1 V þ A s 01 1 p p ¼ A v  q ; (23b) 2 2 p 2 V  A s 02 2 p Mathematical and Computer Modelling of Dynamical Systems 189 A A 1 2 p p 1 2 q q 1 2 p p s t Figure 6. Schematic diagram of the hydraulic actuator. with the bulk modulus of oil β, the piston areas A and A , the offset volumes V and V , 1 2 01 02 and the velocity v ¼ s_ . The volume flows q and q are determined by the valve position p p 1 2 s and the pressure difference, 8 qffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi α A ðÞ s p  p s  0 v1 v s 1 v qffiffi q ¼ (24a) pffiffiffiffiffiffiffiffiffiffiffiffiffiffi α A ðÞ s p  p s < 0; v1 v 1 t v 8 qffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi α A ðÞ s p  p s  0 v2 v 2 t v qffiffi q ¼ (24b) pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi α A ðÞ s p  p s < 0; v2 v s 2 v where α denotes the contraction coefficient, ρ the oil density, p the supply pressure, p the s t tank pressure and A ðÞ s and A ðÞ s are the opening areas of the valve to chamber 1 and v1 v v2 v 2, respectively. The resulting cylinder force F acting on the mechanical part of the system reads as F ¼ p A  p A  F v ; (25) 1 1 2 2 R p where the first part accounts for the pressure force generated by the cylinder and F (v ) R p summarizes the friction forces of the cylinder. It is well known that friction plays an important role in hydraulic actuators. In particular a certain amount of static friction is inevitable for the type of hydraulic cylinders typically used in the considered application. In order to analyse the influence of friction on the dynamic behaviour and the control performance, a static friction model is incorporated into the mathematical model. As can be seen from Figure 7, the friction force is composed of a viscous part r v , a Coulomb V p part r sign(v ) and the static friction r . C p S The static friction phenomenon leads to a switched system behaviour, i.e., if the piston velocity is zero (v = 0) and the absolute value of the sum of the forces on the piston is smaller than r , then the piston sticks and the whole system loses a degree of freedom. S 190 J. Henikl et al. F (v ) R p r v V p −r −r Figure 7. Static friction model. Due to the special choice of the coordinate frame and the basis functions, the detection of zero velocity v ¼ 0 of the ith piston is equivalent to check for φ_ ¼ 0. Furthermore, if p;i i the ith cylinder is sticking, the equations of motion result from Equation (22) by using the constraints φ_ ¼ 0 and φ € ¼ 0. To verify the sticking condition and calculate the absolute i i value of the sum of the forces acting on the piston, the force F of the mechanical part is calculated from the ith row of Equations (20) and (22) N N 1 X X @fðÞ φ € _ _ F ¼ MqðÞ½ i; j q½ j þcqðÞ ; q ½ i þ D½ i; j q½ j þgqðÞ½ i : (26) @φ j¼1 j¼1 Since the functions fðÞ φ are strictly increasing in the operating range, the invertibility of @fðÞ φ =@φ is ensured. Summarizing, the ith piston is sticking if v ¼ 0 and the i p;i i i condition p A  p A  F  r ; (27) 1;i 1;i 2;i 2;i i S;i holds. 4. Parameter identification for a mobile concrete pump In this section, the industrial mobile concrete pump with four joints according to Figure 1 with an operating range of about 40 m is considered. The mathematical model can be derived along the lines of Sections 2 and 3 provided that the following assumptions hold: ● The concrete pump is considered as a planar manipulator. The rotational movement and the lateral beam deflections are neglected. ● The truck’s movement due to the dynamic load of the boom is neglected. ● The boom segments are modelled as Euler–Bernoulli beams. ● The pumping of wet concrete is not considered. However, in view of the last assumption it is worth noting that for the control design the vibrations induced by the pumping are considered as disturbances acting on the system. The model parameters are determined by means of data sheets, design drawings and experimental measurements. In contrast to the parameters of the hydrau- lic components, which are well known from data sheets, the parametrization of the Mathematical and Computer Modelling of Dynamical Systems 191 boom model turns out to be more involved. Besides the geometrical data, the position- dependent mass distribution and the second moment of inertia can be calculated from design drawings. In this context, it is rather difficult to systematically account for externally assembled components on the boom segments like hydraulic cylinders, hydraulic lines and concrete delivery pipes. In particular, their influence on the resulting stiffness of the overall system can hardly be determined analytically. Furthermore, the damping coefficients are unknown and have to be identified by dynamical measurements. 4.1. Experimental set-up In the experimental set-up, in addition to the machine vision system, which will be presented in Section 4.2, the following sensors were used: ● Inclination sensors at both ends of each boom segment. ● Two pressure sensors at each hydraulic cylinder. ● Incremental encoders to measure the cylinder piston movements. Strain gauges attached to the end of each boom segment where the hydraulic actuator is located, in order to determine the beam flection. The strain gauges measure the elongation of the surface of the boom segment. Using the Euler–Bernoulli assumption, the relation between the elongation and the beam flection is given by @ wðÞ x i i ε x ¼ k þ b ; (28) i SG;i SG;i SG;i @x x ¼x i SG;i with the gain k the bias b and the strain gauge position x . In the case of a SG;i SG;i SG;i homogeneous beam, the gain coefficient k is given by the distance of the surface to the SG;i symmetry line of the beam. However, in the case of an inhomogeneous beam it is difficult to determine k analytically. Since in practice it is nearly impossible to attach the strain SG;i gauge to a completely unbent boom segment, the strain gauge signal contains a bias component b , which also depends on the temperature of the beam. Therefore, the gain SG;i and bias coefficients have to be identified by means of measurements. 4.2. Vision system The parameters of the mathematical model of the boom, in particular the mass distribu- tion, the stiffness and damping coefficients and the strain gauge coefficients, could be identified based on static measurements at several configurations of the boom and dynamic experiments with a prescribed excitation of the system. However, due to the large scale of the boom there was the desire to also validate the absolute movement of the boom segments, which cannot be captured by the standard sensors installed. Therefore, a tailored vision system was developed which allows to measure the absolute movement of particular spots on the boom. For this purpose, markers placed at the end of each boom segment and the slewing gear, see Figure 8, were recorded by cameras, and their location in the work space was calculated by means of the respective pixel coordinates on the images. 192 J. Henikl et al. M M M M M 1 2 4 0 3 Figure 8. Markers placed on the boom. In order to capture the large work space of the boom, two 5-mega-pixel greyscale image FireWire cameras were used. The 2452 × 2054 pixel resolution and a frame rate of 15 Hz allow for measurements within a working space of 20 m with an average pixel precision of about 0.01 m/pixel. The open source framework ROS [11] for robotics applications was utilized to record synchronized camera images and to visualize the measurements. After comparing various marker types, including [12–14], circular markers without encoded identifiers were selected. The detection of the circles and their projected elliptical shapes was done by using the standard computer vision library OpenCV. With the algorithm presented in [15] and [12], the precision of the marker position estimation can be increased up to sub-pixel level. This enables the detection of boom movements of even 2 mm. In addition to the markers on the boom, further markers located at fixed positions were used for offline estimation of static camera poses. In order to start the tracking algorithm, all measurement points must be initialized in the camera image. 4.2.1. Feature detection In order to speed up the marker detection process, the algorithm has an initial feature location and predicts the next location throughout the process. For this reason, the marker detection has to be applied only to a small window around the predicted feature location. During the initialization phase, pixel locations and the estimated size of the circular feature have to be defined. The detection algorithm then detects multiple ellipse hypoth- eses in each predicted search window. For this purpose, a function which uses the estimated location and the expected radius sorts all detected ellipses within the search window and places the most likely one first. This ellipse is then used as the basis for the next search window. Figure 9 shows multiple search windows with the most likely ellipse in red. 4.2.2. 3D data generation The geometric relation between the marker position in the workspace and their projected image coordinates are given in this section. The transformation of a point p ¼½p ; p ; p  , defined in the workspace coordinate frame 0 x y z , to the W;x W;y W ;z W W W W W Mathematical and Computer Modelling of Dynamical Systems 193 M M M M M M 2 3 4 0 1 2 (a) (b) Figure 9. Search windows with detected edges in green, ellipses in yellow and the most likely marker ellipse in red. (a) Left camera. (b) Right camera. respective coordinates p ¼½p ; p ; p  in the camera coordinate frame 0 x y z , C;x C;y C;z C C C C 31 which is defined by the camera centre, is given by the translation vector d 2 R and 33 the rotation matrix R 2 R , W W p ¼ R p þ d : (29) C C W C In the literature, this transformation is typically described by the so-called external camera matrix, see e.g., [16]. Assuming a pinhole camera model, the projection of the point p onto the 2D image frame 0 u v is given by I I I "# C;x f þ o p ; u u I u C;z p ¼ ¼ ; (30) I C;y p ; f þ o I v v v C;z with the internal camera parameters f ¼ f =s ; f ¼ f =s ; o and o . This projection is u u v v u v typically referred to as intrinsic camera matrix, see e.g., [16]. The focal length f describes the distance of the camera centre to the sensor chip or the image frame, respectively, s and s represent the effective physical dimensions of one pixel on the camera chip and the distance between the image centre and the origin of the image frame is given by o and o . u v In order to compensate for lens distortions and to be able to use a simple pinhole camera model, every frame captured has to be preprocessed first to get an undistorted image before it is further used in any detection process, see e.g., [16]. The parameters for the compensation of the lens distortions as well as the internal camera parameters were estimated by means of an offline calibration procedure. The calibration was done using a checker-board pattern of known size and the ROS camera calibration node. For the estimation of the three linearly independent quantities parame- W W trizing the rotation matrix R and the three coefficients of the translation vector d ,an C C iterative method based on the Levenberg–Marquardt [17] optimization was utilized, which is implemented in the before mentioned OpenCV library. For this purpose, the workspace coordinates of at least three known reference points p and their visual appearance in the image frame coordinates p are required. Thus, with the transformation (29) and the projection (30), six non-linear equations are given to determine the six unknown coeffi- W W cients for R and d . With the use of additional reference points, the error due to the C C measurement uncertainty of the reference coordinates could be minimized. In Figure 10, five reference points p for j ¼ 1; ... ; 5, used for the estimation of the camera pose, are W ;j shown. The coordinates were measured manually by a laser distance sensor. The reference 194 J. Henikl et al. W,3 W,2 e e p p W,4 W,5 W,1 Figure 10. Reference points. e e points p for l ¼f1; 2; 3; 4g and p for r ¼f1; 2; 3; 5g were used for the left and the W ;l W ;r right camera, respectively. With these results, the relation of the coordinates of a marker centre p described in the workspace coordinate frame and the respective coordinates p is given. The goal here is, however, the inverse problem. A 3D point linked to a 2D image point has to be found. This can be done by creating a straight line, see e.g., [18] that intersects the camera centre and the image frame at the respective 2D image point. Considering Equation (29) and p ¼ 0, the camera centre described in the workspace frame is given by c W W p ¼ R d : (31) W C C With Equation (30) and p ¼ f , the 2D image point at the image frame in front of the C;z m m m camera, defined by the pixel coordinates of the marker centre p ¼½p ; p  ,is I I ;u I ;v located at 2 3 1 m p  o f I;u 6 7 p ¼ 1 m (32) 4 5 p  o f f I;v and m W 1 m W 1 W ^ ¼ðR Þ p ^ ðR Þ d ; (33) W C C C C respectively. With this and Equation (31), the direction of the line is given by 2 3 1 m p  o I;u 6 7 m c W 6 7 v ¼ p  p ¼ f R 1 m ¼ f v (34) W W W W C 4 5 p  o f I;v and thus, the line could be described by the equation Mathematical and Computer Modelling of Dynamical Systems 195 l c p ðλÞ¼ p þ λv (35) W W for λ 2 R. The final 3D marker location p can be computed under the assumption that the boom and the marker moves only in the plane defined by n p ¼ λ ; (36) with the normal vector n to the plane and a scalar λ defining the position of the plane. The normal vector n and the scalar λ were determined by means of the reference points lying in the plane, see Figure 10. With this, the scalar λ has to be found that defines the c m intersection of the line (35) with the plane (36), p þ λ v ¼ p . Considering τ W W W T c n ðp þ λ v Þ¼ λ ; λ is given by τ W  τ T c λ  n p λ ¼ ; (37) n v and therefore p reads as T c T c λ  n p λ  n p m c  W c  W p ¼ p þ v ¼ p þ v : (38) W W W W W T T n v n v W W 4.3. Identification procedure Due to the model complexity and the high number of system parameters, the identification has to be performed in several steps. The following procedure was applied to the simulation model of the mobile concrete pump. First, a fine tuning of the mass distribution was obtained by means of static experi- ments. For this purpose, the cylinder forces respectively the chamber pressures and the inclination at the beginning of each boom segment were measured at different angles φ ¼ φ  of the first joint and a sprawled boom, i.e., φ ¼ φ ¼ φ ¼ 0. For the pre- 1;m 1 2 3 4 scribed constant angles φ  , the corresponding measured cylinder forces and inclinations 1;m are referred to as F ðφ  Þ and ψ  φ  ; i ¼ 1; ... ; N. In order to compare the i;m 1;m b;i;m 1;m measurement results with the simulation model, the stationary solution of Equations (20) and (22) is determined, i.e., gðq Þ¼ Q and, given φ ¼ φ ¼ φ ¼ 0, solved for 2 3 4 b b o o h ; ... ; h ; h ; ... ; h and F ; ... ; F as a function of φ  . In general, the inclination at the 1 N 1 1 N 1 N beginning of the boom segment reads as i i1 X X ψ ¼ φ þ w ðL Þ: (39) n n b;i n n¼1 n¼1 Thus, for φ ¼ φ ¼ φ ¼ 0 and w (x ) according to Equation (17), the stationary 2 3 4 n n relation of Equation (39) takes the form i1 b b o o ψ ¼ φ þ h η ðÞ L þ h η ðÞ L : (40) b;i 1 n n n n n n n¼1 196 J. Henikl et al. With Equation (40), the stationary cylinder forces F ; i ¼ 1; ... ; N, can be parame- trized in terms of the respective inclination, i.e., F ðψ  Þ. Now the nominal mass distribu- i b;i tion of the boom segments σ ðx Þ; i ¼ 1; ... ; N is scaled with a multiplicative constant i i ; i ¼ 1; ... ; N, which serves as a parameter to be identified. For this, the least squares problem min F ðψ  Þ F ðψ  Þ (41) i i;m b;i;m b;i;m has to be recursively solved beginning at the last boom segment i = N. In order to identify the stiffness and the damping behaviour of the boom, dynamic experiments with a defined excitation of the system were performed in the next step. With the measurements of the cylinder forces, the strain gauges and the marker movements of the machine vision system, the mean equivalent Young’s modulus E as well as the damping coefficients d and d were determined such that the simulation model matches φ w the occurring oscillation frequencies and the decay of the envelopes. The strain gauge coefficients k and b were also identified by means of static measurements, where a SG;i SG;i similar procedure as for the adjustment of the mass distribution was employed. 5. Experimental results The following subsection presents the validation of the simulation model by means of experimental results. The numerically calculated basis functions (14) and (15), employed in the simulation model, are illustrated in Figure 11.In Figure 12, the results of static measurements in comparison to the simulation model are given. The cylinder forces and beam flections have been measured at different angles φ of the first joint and a sprawled boom with joint angles φ , φ and φ close to 0. The cylinder forces are normalized to the 2 3 4 1 1 0 0 x = x 1 m −1 −1 0 0.2 0.4 0.6 0.8 1 00.2 0.4 0.6 0.8 1 x /L x /L 1 1 2 2 1 1 0 0 −1 −1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x /L x /L 3 3 4 4 1st basis function 2nd basis function Figure 11. Basis functions of the boom segments normalized to the beam length and to the maximum function value. η (x )/max|η (x )| η (x )/max|η (x )| 3 3 3 3 1 1 1 1 η (x )/max|η (x )| η (x )/max|η (x )| 4 4 4 4 2 2 2 2 Mathematical and Computer Modelling of Dynamical Systems 197 0. 6 0.5 0. 4 0.3 0. 2 0.1 0. 5 0.3 0. 3 0.2 0.1 0. 1 500 400 250 150 020 40 60 80 020 40 60 80 ϕ [°] ϕ [°] 1 1 Simulation Measurement Figure 12. Static measurements at a sprawled boom: cylinder forces normalized to the maximum forces and strain gauge measurements as a function of the first joint angle. maximum forces, which are determined by the supply pressure and the piston area, F ¼ p A . i;N s 1:i In order to validate the dynamic system behaviour, the system was excited by a sudden release of a 90 kg load at the end of the boom. In Figure 13, the comparison of the measured and simulated response of the normalized cylinder forces, the strain gauge signals and relative movements of the cylinder pistons at a typical boom configuration on the construction site is shown, see Figure 1. The stationary deviations of the strain gauge signals can be explained by thermal expansion due to different experiment days. Finally, in Figure 14 the measurement results of the machine vision system are presented. The simulated and tracked movement of the markers, described in the inertial frame, is illustrated. It can be seen that the amplitudes of the oscillations agree very well. With this, an essential requirement for the mathematical model, which serves as a basis for the design of active damping control strategies, is fulfilled. Looking at the absolute deviations to the inertial frame, one has to bear in mind that already small deviations of the angle w ðL Þ at the end of the first boom segment have a big influence on the absolute marker positions. Furthermore, due to lens distortions at the outer image section as well as the restricted validity of the assumption that all markers move in the same plane, the precision of the absolute machine vision measurements is limited. However, considering the size of the boom, the biggest vertical and horizontal deviations of about 60 cm for y , m 4 and 15 cm for x can be considered small in relation to the horizontal range of 27 m. m;3 ε [μm/m] ε [μm/m] F /F F /F 3 3 3,N 1 1,N ε [μm/m] ε [μm/m] F /F F /F 4 2 4 4,N 2 2,N 198 J. Henikl et al. 0.62 0.56 0.58 0.52 0.48 0.54 0.4 0.44 0.4 0.3 0.36 0.2 –50 0.6 0.3 0.4 0.2 0.1 0.2 0.6 0.8 0.6 0.4 0.4 0.2 0.2 05 10 15 20 0 510 15 20 t [s] t [s] Simulation Measurement Figure 13. Load release at the end-effector: cylinder forces, strain gauge signals and cylinder piston movements. Since in view of the envisaged controller design the main goal of the model is to capture the essential dynamic behaviour of the mobile concrete pump, the impreciseness in terms of the absolute deviations is of minor importance. Figure 14 shows the relatively high resolution and low noise level of the signals, taking into account the large work space. Figure 15 depicts the relative vertical movement of a marker attached to the slewing gear. It can be seen that the vision system is even able to capture the 2 mm amplitude of the reaction of the truck caused by the sudden release of the load. With the use of the effective ellipse search algorithm, the image processing time for each frame is far lower than the sampling time defined by the frame rate. Hence, it would even be possible to use this system for real-time applications. Δs [mm] Δs [mm] ε [μm/m] ε [μm/m] F /F F /F p,3 p,1 3 1 3 3,N 1 1,N Δs [mm] Δs [mm] ε [μm/m] ε [μm/m] F /F F /F p,4 p,2 4 2 4 4,N 2 2,N Mathematical and Computer Modelling of Dynamical Systems 199 14.7 6.77 6.75 14.5 6.73 14.3 27.5 21.6 21.55 27.3 21.5 21.45 27.1 4.1 5.3 5.2 4.05 5.1 2.6 −2.4 2.4 −2.8 2.2 −3.2 05 10 15 20 05 10 15 20 t [s] t [s] Simulation Left camera Right camera Figure 14. Load release at the end-effector: horizontal and vertical distance of the markers to the inertial frame. −1 −2 −3 02468 10 12 14 16 18 20 t [s] Figure 15. Load release at the end-effector: relative vertical movement of a marker attached to the truck. 6. Conclusions In this paper, a systematic approach for the mathematical modelling of large-scale light- weight manipulators with hydraulic actuation was proposed. It was shown that a suitable choice of the basis functions for the approximation of the elastic bending of the beams has significant advantages for the incorporation of static friction into the model. Furthermore, simulation and measurement results of a mobile concrete pump y [m] y [m] x [m] x [m] m,3 m,1 m,3 m,1 Δ [mm] ym,0 y [m] y [m] x [m] x [m] m,4 m,2 m,4 m,2 200 J. Henikl et al. show the accuracy and the feasibility of the proposed modelling approach already for a small number of basis functions. A procedure for the identification of the uncertain model parameters is given. For this purpose, a tailored machine vision system was presented for the measurement of the absolute boom movement. It was demonstrated that this system enables the tracking of motions with amplitudes of approximately 2 mm even for a large observation area. Acknowledgements The authors want to express their gratitude to Helmut Brunner for the implementation of the vision system algorithm as well as the recording and evaluation of the camera images. Funding Part of the research of this work was supported by the company TTControl and the Austrian Research Promotion Agency (FFG) within the project MOBILE [Project no. 2426433]. Supplemental data Supplemental data for this article can be accessed here. Notes 1. VICON: http://www.vicon.com 2. For the sake of readability, the time dependence of the degrees of freedom is not explicitly stated in this paper. 3. Allide Vision: Pike F-505, http://www.alliedvisiontec.com 4. ROS (Robot Operating System): http://www.ros.org 5. OpenCV (Open Source Computer Vision): http://opencv.willowgarage.com/wiki/ 6. ROS Camera Calibration: www.ros.org/wiki/camera_calibration 6 1 7. With the linear thermal expansion coefficient for steel of α  11 ... 13  10 K a T ;Steel temperature difference over night of 10 K yields an expansion of about 110 ... 130 μm=m: 8. A video illustrating the capturing of the vertical motion of the markers is available on http:// www.acin.tuwien.ac.at/fileadmin/cds/videos/Motion_Capturing_Mobile_Concrete_Pump.mp4 References [1] A.A. Shabana, Dynamics of Multibody Systems, Cambridge University Press, New York, 2005. [2] H. Bremer, Elastic Multibody Dynamics: A Direct Ritz Approach, Springer, Heidelberg, 2008. [3] M. Schneider, Modellbildung, Simulation und nichtlineare Regelung elastischer, hydraulisch angetriebener Großraummanipulatoren, VDI Fortschrittsberichte Reihe 8, VDI Verlag, Düsseldorf, 1999. [4] H. Aschemann, O. Sawodny, A. Bulach, and H.P. Hofer, Model based trajectory control of a flexible turntable ladder, Proceedings of the American Control Conference Vol. 2, Anchorage, AK, 8–10 May 2002, pp. 921–926. [5] S. Lambeck, O. Sawodny, and E. Arnold, Trajectory tracking control for a new generation of fire rescue turntable ladders, Proceedings of the 2nd IEEE Conference on Robotics, Automation and Mechatronics, Bangkok, 1–3 June 2006, pp. 847–852. [6] M. Dibold, J. Gerstmayr, R. Stadlmayr, H. Irschik, and K. Schlacher, Dynamics of multibody systems including hydraulic actuators and feedback control, Proceedings of the 3rd European Conference on Structural Control Vol. 2, Vienna, 12–15 July 2004, pp. 39–42. [7] A. De Luca and B. Siciliano, Closed-form dynamic model of planar multilink lightweight robots, IEEE T Syst. Man Cy. 21 (1991), pp. 826–839. doi:10.1109/21.108300. [8] J. Henikl, W. Kemmetmüller, and A. Kugi, Modeling and control of a mobile concrete pump, Proceedings of the 6th IFAC Symposium of Mechatronic Systems, Hangzhou, 10–12 April 2013, pp. 91–98. Mathematical and Computer Modelling of Dynamical Systems 201 [9] R. Mahony, V. Kumar, and P. Corke, Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor, IEEE Robot. Autom. Mag. 19 (2012), pp. 20–32. doi:10.1109/ MRA.2012.2206474. [10] J. Henikl, W. Kemmetmuüller, and A. Kugi, Modeling and simulation of large-scale manip- ulators with hydraulic actuation, Proceedings of the 7th Vienna Conference on Mathematical Modeling, Vienna, 14–17 February 2012, pp. 780–785. [11] M. Quigley, K. Conley, B. Gerkey, J. Faust, T. Foote, J. Leibs, R. Wheeler, and A.Y. Ng, ROS: An open-source robot operating system, ICRA Workshop on Open Source Software, Kobe, 12–17 May 2009. [12] X. Zhang, S. Fronz, and N. Navab, Visual marker detection and decoding in AR systems: A comparative study, Proceedings of the 1st International Symposium on Mixed and Augmented Reality, Darmstadt, 30 September–1 October 2002, pp. 97–106. [13] L. Naimark and E. Foxlin, Circular data matrix fiducial system and robust image processing for a wearable vision-inertial self-tracker, Proceedings of the 1st International Symposium on Mixed and Augmented Reality, Darmstadt, 30 September–1 October 2002, pp. 27–36. [14] F. Ababsa and M. Mallem, A robust circular fiducial detection technique and real-time 3D camera tracking, J. Multimedia. 3 (2008), pp. 34–41. doi:10.4304/jmm.3.4.34-41. [15] M. Zillich and J. Matas, Ellipse detection using efficient grouping of arc segments, Proceedings of the 27th Workshop of the Austrian Association for Pattern Recognition (OAGM/AAPR), Laxenburg, 5–6 June 2003, pp. 143–148. [16] E. Trucco and A. Verri, Introductory Techniques for 3-D Computer Vision, Prentice Hall PTR, Upper Saddle River, NJ, 1998. [17] D. Grest, T. Petersen, and V. Krüger, A comparison of iterative 2D-3D pose estimation methods for real-time applications,in Image Analysis, A.B. Salberg, J. Hardeberg, and R. Jenssen, eds., Springer, Heidelberg, 2009, pp. 706–715. [18] C. Wohler, 3D Computer Vision: Efficient Methods and Applications, Springer, Heidelberg,

Journal

Mathematical and Computer Modelling of Dynamical Systems – Taylor & Francis

Published: Mar 4, 2015

Keywords: flexible multibody system; hydraulic cylinder; parameter identification; machine vision

There are no references for this article.