Analytical methodology for the analysis of vibration for unconstrained discrete systems and applications to impedance control of redundant robots

This paper presents a general methodology for the analysis and synthesis of a positive semi-definite system described by mass, damping and stiffness matrices that is often encountered in impedance control in robotics research. This general methodology utilizes the fundamental kinematic concept of rigid-body and non-rigid-body motions of which all motions consist. The rigid-body mode results in no net change in the potential energy from the stiffness matrix of the multiple degree-of-freedom (DoF) discrete mechanical system. Example of an unconstrained discrete mechanical system is presented to illustrate the theoretical principle as applied in obtaining the free and forced vibration responses, as well as the dynamic characteristics of the system in natural frequency, ωn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega_n$$\end{document} and damping ratio, ζ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta$$\end{document}. In addition, the methodology is applied to the impedance control of redundant robots. The rigid-body mode is equivalent to the motions of a redundant robot which result in no net change in potential energy, also called the zero-potential or ZP mode, of impedance control. Example of a redundant robot is used to demonstrate the application of the methodology in robotics. The dynamic characteristics of ωn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _n$$\end{document} and ζ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta$$\end{document} in the modal space are analyzed, which can be synthesized to modulate the damping of the system analytically.


Introduction
When impedance control is utilized in robotics research, it is often realized by researchers in the field that the elements of the matrix, especially the damping matrix to impart proper damping to achieve desired dynamic response, are obtained using trials and error without available analytical tools. This is largely due to the multiple degrees of freedom (DoF) of the robotic system, as manifested by coupled stiffness and damping matrices that make the analysis and synthesis challenging. We present in this paper a new and general methodology for the analysis and synthesis of a positive semi-definite robotic system described by mass, damping and stiffness matrices in impedance control. This general methodology grew out of a novel research result in vibration analysis that involves general mass, damping, and stiffness matrices of unconstrained mechanical systems.
The analysis of the vibration response of discrete mechanical systems with masses, dampers, and springs has been developed for constrained discrete systems [1]. However, a systematic methodology for the solution of unconstrained systems with multiple degrees of freedom (DoF) has not yet been presented. In this paper, a new and general methodology is presented for unconstrained discrete mechanical systems of multiple DoF by utilizing the kinematic concept of rigid-body (RB) and non-rigidbody (NRB) motions. The rigid-body motions refer to the motion of such multiple DoF system with no net change in potential energy from the interconnected springs; that is, the rigid-body motions belong to the null space of the stiffness matrix. The unconstrained systems possess rigid-body mode(s) that must be removed before solving the vibration problem and determining the dynamic response.
The methodology is also applied to the analysis of impedance control of redundant robots, in which the redundancy resembles the unconstrained rigid-body mode(s). Such rigid-body mode(s), or better described as the zero-potential-energy mode(s), results in no net change in the potential energy owing to the stiffness in the impedance control of redundant robotic manipulators. This method enables the analysis of dynamic response of impedance control of redundant robots with given mass, damping and stiffness matrices. Furthermore, it enables one to choose the elements of stiffness and damping matrices to modulate the dynamic responses of the robot with proper damping ratios by employing this analytical method.
The main contribution of this paper is to provide a systematic methodology, based on the vibration theory and kinematics of motion, to determine the dynamic characteristics of a system under impedance control. Furthermore, one can use this methodology to analytically modulate and optimize the damping characteristics of the impedance control by gaining physical insights into the dynamic characteristics of the mechanical or robot system using the analytical methodology.
This paper can logically be separated into two parts. Part I presents the theory of analysis and methodology using the eigenvalue analysis. Part II presents the application of the theory in robotic impedance control on a redundant robotic manipulator.

Related work
There has been work by many authors on vibration analysis of mechanical systems that are constrained to an inertial frame of reference, such as in [1], among others, where the systems are conservative and positive definite. For the case of constrained non-conservative systems, such as positive-definite mass-damper-spring systems of multiple degrees of freedom, the linear system theory [2] can be directly applied to solve them to obtain their dynamic response. If such systems are unconstrained (positive semi-definite), they are no longer readily solvable, unless the methodology presented in this paper is employed.
Impedance control [3] of robotic manipulators deals with the interaction of the robot with its environment, and is modeled as a multi-dimensional mass-damperspring system, similar to the discrete systems in mechanical vibration theory. The methodology presented here becomes particularly useful for the case of redundant manipulators, which require a special treatment, as in non-conservative unconstrained discrete mechanical systems. This methodology allows us to obtain the dynamic response of the robot, which also lets us analyze the effect of different damping parameters of the system [4] once it has been appropriately mapped into the joint space [5,6] (as most of robotic tasks are defined in Cartesian space [7]) to damp out undesired joint vibrations in a theoretically sound manner. Research results on selecting damping parameters for impedance control were presented, for example in [8], but assumed a special case of diagonalization of the system with both damping and stiffness, which in general is not the case.

Motivation and essence of the method
Impedance control of systems with multiple degrees of freedom (DoF), represented by Eq. (1), requires knowledge of the dynamic characteristics of the system for good control. 1 Equation (1) represents a multiple-DoF robotic system under impedance control with the {M, C, K} matrices. The challenge of such a system is the non-diagonal (or coupled) {M, C, K} matrices that make the system equation difficult to decouple to obtain the dynamic characteristics of a scalar system, like the one-DoF system of mẍ + cẋ + kx = f . (Even when one chooses diagonal C and K matrices, the M matrix is still not diagonal for any robotic system.) Furthermore, the elements in damping C are often chosen by trial-anderror or from experience.
Take the case of impedance control of a Franka Emika Panda robot at a given configuration with a prescribed K = diag ([2000, 600, 2000, 100, 100, 100]) in Cartesian space, one can choose a diagonal damping matrix of C 1 = diag([5, 3, 5, 1, 1, 1]) (in SI units). What often follows is a trial-and-error approach to pick the elements of C to achieve a good dynamic characteristics by running the robot and measuring its impulse or step response. An experimental response of this set of {M, C 1 , K} (with the mass matrix M given by the robot system library) is plotted in Fig. 1, with large under-and over-shoot. By employing the methodology in this paper to obtain dynamic responses and understand the dynamic characteristics in the modal space in which the system Eq. (1) can be decoupled to 6 independent single DoF, we can choose the following C 2 = diag ([22.5, 23.5, 45, 2, 12, 1]) . This can be easily done by following the proposed method using a range of values of individual elements in C to study their trend of influence on the damping ratios using automated code. The much improved dynamic response of experimental results on Panda robot at the same configuration is plotted in Fig. 1 for comparison.
It follows from the theoretical analysis in this paper, the damping ratios and natural frequencies of the decoupled system in the modal space are shown in Table 1. The system with C 2 produces larger damping ratios; whereas, C 1 produces 3 small damping ratios. This explains the dynamic response of low damping using C 1 , and improved response using C 2 in the experimental results of responses in Fig. 1. It is noted that the response in the Cartesian space is a linear combination of each response in the modal space, according to the well-known "expansion theory" in vibration analysis [1].
The key is to choose the corresponding elements of the damping matrix using the analytical methodology to understand the fundamental dynamic behavior and improve the dynamic response of the system. This is accomplished by the novel methodology based on vibration theory and kinematics of motion, presented in "Theoretical background" section. Examples of a mechanical system and a redundant robot system are presented with illustration in the two sections following the theory.

Theoretical background
An unconstrained discrete mechanical system has freedom to move without being constrained to an inertial frame, such as the constraint of a wall or a rigid structure. Such unconstrained system has rigid-body motion, as if the entire system were a rigid-body moving without relative motions between the discrete elements of the system [9]. This rigid-body mode of motion is dictated by rigidbody mechanics, with a frequency of oscillation of zero. We will denote the rigid-body mode of motion of an n-DoF system as u 0 , which is a n × 1 vector of constant elements, with the associated frequency of ω 0 = 0.
First of all, let us consider the general equation of motion of a n-DoF system in Eq. (1) where M is the mass matrix, C is the damping matrix, K is the stiffness matrix of the n-DoF dynamic system described by the n × 1 state vector q(t) = [q 1 (t) q 2 (t) · · · q n (t)] T subject to an external force Q . The initial displacements and velocities are given as q 0 = q(0) and q 0 =q(0) , respectively.
The characteristic of an unconstrained, positive semidefinite system is a singular stiffness matrix, K . The system in Eq. (1) becomes a positive semi-definite system. In addition, the rigid-body mode u 0 belongs to the null space of K [10]; that is, As an example, the rigid-body mode of the unconstrained system shown in Fig. 2 is u 0 = [1 1 1] T when all three masses of the three-DoF mechanical system move synchronously, as if they were a rigid body. This concept of "rigid-body mode" also applies to a redundant system, such as a redundant robot in which the number of DoF (joints) is larger than the intended DoF of the system, Experimental results of impedance control on a Panda robot. The responses are performed with an initial displacement in the y direction, released with zero initial speed, corresponding the damping matrices of C 1 and C 2 , respectively typically 6 for a spatial robot with 6 DoF in the Cartesian space. Therefore, Eq. (2) is extended to include such redundant robot system by defining the "zero-potentialenergy" criterion as For the sake of presentation in this paper, however, the notation of rigid-body mode is used, although it is interchangeable with the "zero-potential-energy" criterion throughout this paper.
This rigid-body mode defined in Eq.
(2) can be normalized with respect to the mass matrix by u T 0 M u 0 such that The corresponding non-oscillatory frequency of this rigid-body mode is In the following, we provide the procedures to solve the unconstrained or redundant systems with: (I) remove the redundancy, or the unconstrained rigid-body mode; (II) formulate the positive-definite system; (III) determine the response of the system; and (IV) obtain the overall response of vibration by combining the RB and NRB solutions. An example is illustrated with an unconstrained redundant mechanical system of three degrees of freedom. Finally, the application of the theory to impedance control of redundant robots is presented with an example.

Remove the redundancy
First of all, we must remove the redundancy, or the unconstrained rigid-body motion, of the positive semidefinite system described in Eq. (1).
By the expansion theorem, a general response of vibration can be expressed as a linear combination of all the bases (modes, or modal vectors) in the following where q RB (t) and q NRB (t) denote the rigid-body and non-rigid-body components of the state vector q(t) . All modal vectors, u i , including u 0 , are normalized such that u T i Mu i = 1 for i = 0, 1, 2, . . . , (n − 1) . By the orthogonality property, we have Next, in order to determine the oscillatory motions (or the non-rigid-body, NRB, motion), we remove the known rigid-body mode by imposing upon the state vector q a constraint matrix, S , to ensure that all solutions of q(t) are free of the rigid-body mode, based on Eq. (6). By the expansion theorem, any response q(t) is a linear combination of all modal vectors, including the rigid-body mode. If a response q(t) is free of the rigid-body mode, it will be orthogonal to the rigid-body mode, u 0 ; that is, Equation (9) can be used to equate a chosen q i as a function of the rest of (n − 1) elements in q by the following equation Therefore, the constraint matrix, S , can be defined and expressed as q = S q ′ where where the constraint matrix and the reduced state vector q ′ are (6) s i q i = s 1 q 1 + s 2 q 2 + · · · + s n q n = 0 (10) The new state vector q ′ = [q 1 q 2 · · · q i−1 q i+1 · · · q n ] T is the reduced state vector of q by removing q i . For example, if we choose to eliminate q n , Eq. (10) will become The constraint matrix in Eq. (11) will become q = S q ′ or Employing the mapping of q = S q ′ in Eqs. (11) or (13), Eq. (1) can be pre-multiplied by S T to become where the matrices in the reduced (n − 1) space with the new state vector q ′ are defined as

Formulate the positive-definite system
Once the rigid-body mode has been removed, the equation of motion of the reduced (n − 1)-DoF, as derived above, is formulated as where q ′ (t) = [q 1 q 2 · · · q i−1 q i+1 · · · q n ] T is the reduced vector with (n − 1) state variables, defined by q = Sq ′ in (12) q n = − s 1 s n q 1 + s 2 s n q 2 + · · · + s n−1 s n q n−1 Eqs. (11) or (13). The equation of motion in (14) with the M ′ , C ′ and K ′ matrices is now a positive definite system. 2 The solution of Eq. (14) requires the excitation Q and the initial conditions q ′ (0) and q ′ (0) , which have to be derived from the given initial conditions q(0) and q(0) in the physical space. Based on Eq. (6), the initial condition q(0) can be broken into RB (rigid-body) and NRB (nonrigid-body) components where α 0 is a constant. Pre-multiply Eq. (15) by u T 0 M to obtain Note that u T 0 Mu 0 = 1 and u T 0 Mq NRB (0) = 0 based on the orthogonality in Eq. (7). The RB part of the initial displacement is Likewise, the RB part of the initial speed can be expressed as Thus, the NRB part of the initial displacement and speed for the equation of motion (14) in the reduced (n − 1) space, free of the rigid-body mode, are The initial conditions of the reduced state vector q ′ (0) and q ′ (0) will be the corresponding elements in q NRB (0) and q NRB (0) from Eqs. (18) and (19), respectively. For example, if we choose to eliminate q n , the first (n − 1) elements of q NRB (0) and q NRB (0) in Eqs. (18) and (19) should be used as the initial conditions of the reduced state vectors q ′ (0) and q ′ (0) , respectively.

Determine the response of the system
Finally, the solution of vibration analysis of an unconstrained system consists of both (15) 1. Rigid-body mode: Apply Eqs.
(2) and (4) to find the normalized rigid-body mode, u 0 , corresponding to ω 0 = 0. 2. Non-rigid-body mode: Remove the rigid-body mode from the solution of q(t) by using the constraint matrix S to obtain the reduced equation of motion for (n − 1) DoF in Eq. (14). Solve this equation with the given excitation ( S T Q , for forced vibration response) and the initial conditions of the reduced vectors, q ′ (0) and q ′ (0) , from the NRB part of the initial conditions q NRB (0) and q NRB (0) in Eqs. (18) and (19), respectively.
In the following subsections, the two parts of the solution for free and forced vibration responses will be discussed.

Determine the rigid-body motion (RB)
Equations (2) and (4) render the rigid-body mode which is normalized, called u 0 . Since this RB motion depends on K , we need to determine it using the dynamics of the system. With the normalized rigid-body mode, we can write the dynamic response of the rigid-body as follows When u 0 is normalized, as in Eq. (4), Eq. (21) can be simplified to Equation (22) is a scalar equation if there is only one rigid-body mode. 3 Next, we formulate the equations to determine the rigid-body motion with free or forced vibration.

RB motion: free vibration analysis
Free vibration is a response to only the initial conditions. When Q = 0 , the differential equation of motion for β(t) in Eq. (22) becomes The solution can be obtained from the following equation with the initial conditions: The solution of response to Eq. (23) is where a = u T 0 Cu 0 � = 0 is a positive scalar for a positive definite C matrix.
If the system is without damping or a = u T 0 Cu 0 = 0 , we can re-write Eq. (23) as follows The solution of response to Eq. (25) is Once the response β(t) is found, the RB part of the free vibration response can be obtained by Eq. (20), q RB (t) = u 0 β(t) . Note Eqs. (24) and (26) can result in net steady-state RB motion of the system as a whole, as a consequence of the RB part of the initial conditions. (22) can be employed to find the forced vibration response. The solution of this one-DoF system from Eq. (22) can be obtained using the inverse Laplace transform Once the response β(t) is found using Eq. (27), the RB part of the forced vibration response can be obtained by Eq. (20), q RB (t) = u 0 β(t).

Determine the oscillatory motions (NRB)
Once the rigid-body mode has been removed, the equation of motion of the reduced (n − 1)-DoF positivedefinite system in Eq. (14) can be solved to obtain the oscillatory motions (NRB). Because q ′ (t) is free of RB motion, the initial conditions employed for the solution in Eq. (14) must also be free of the RB motion.
Note that the reduced state variables in q ′ (t) still represent the physical generalized coordinates. The differences between q(t) and q ′ (t) include: • The reduced state q ′ (t) has (n − 1) independent generalized coordinates, with the RB motion removed from the n-DoF q(t) , by imposing the constraint matrix S.
• The reduced state q ′ (t) should be free of the RB motion. This has a major implication on how to define the initial conditions when solving for the NRB motion using Eq. (14) with the given initial conditions q(0) and q(0) in the original n-DoF space.
Note that if this system has no damping; i.e., the damping matrix C in Eq. (1) does not exist, a modal matrix in the reduced (n − 1) × 1 space can be found, by the modal analysis, and mapped back to the n-DoF system to find a full n × n modal matrix [1] with the addition of the rigidbody mode, u 0 . This can be treated as a special case of the general treatment which is presented in the following. Next, we formulate the equations to determine the non-rigid-body (NRB) motion with free or forced vibration.

NRB motion: free vibration analysis
The initial conditions are needed to find the free vibration response of NRB motion in the q ′ space. The general RB motion is defined in Eq. (20). The RB part of the initial conditions are in Eqs. (16) and (17), with the RB part of free vibration response derived in "Determine the rigid-body motion (RB)" section. The NRB part of the initial conditions are derived in Eqs. (18) and (19). The corresponding (n − 1) elements of these initial conditions in the reduced q ′ space are used to solve for the free vibration response q ′ (t) based on the equation of motion in (14).
The following two types of systems are discussed: (i) conservative systems without damping, and (ii) nonconservative systems with damping.
i. Conservative systems without damping: Apply the modal analysis for the conservative system [1] to the reduced systems in Eqs. (14), without the damping matrix C , to obtain the (n − 1) × (n − 1) modal matrix, U ′ . Use the corresponding elements of q ′ and q ′ in q NRB (0) and q NRB (0) to form the initial conditions in the reduced space q ′ (0) and q ′ (0) . The initial conditions in the modal space, each a (n − 1) × 1 vector, will be: The solution of the modal coordinates are thus where r = 1, 2, · · · , (n − 1) . Once the vector of the modal coordinates η(t) = [η 1 (t), η 2 (t), · · · , η n−1 (t)] T is obtained, the (n − 1) × 1 reduced state vector is Finally, the NRB part of the free vibration response is obtained from the mapping using the constraint matrix in equation (11) ii. Nonconservative systems with damping: The nonconservative system in Eq. (14) has (n − 1) × (n − 1) mass, damping, and stiffness matrices M ′ , C ′ , and K ′ . Equation (14) can be rearranged in the form of the linear system equation in the following

and the A and B matrices are
In order to find the solution to the system described by Eq. (29), we employ the dual eigenvalue analysis for the nonconservative system. The solution is where and Y are the right and left eigenvectors, respectively, and Y T X = I . Thus, we can substitute the results into Eq. (30) to obtain The free vibration response to the initial conditions without excitation can be obtained from Eq. (31) directly. The response is The eigenvectors X and Y are obtained from the dual eigenvalue analysis together with the eigenvalues in the matrix. The initial conditions are given with x(0) = [q ′ (0)q ′ (0)] T , using the corresponding elements of q ′ and q ′ in q NRB (0) and q NRB (0) to form the initial conditions in x(0) . The NRB part of the free vibration responses, including generalized coordinates and generalized velocities, are found in Finally, the NRB part of the free vibration response is obtained from the mapping using the constraint matrix in Eq. (11) NRB motion: forced vibration analysis Equation (14) is an (n − 1)-DoF positive-definite system. The solution of forced vibration response can be obtained by employing the standard procedures of the modal analysis if the damping C does not exist. For nonconservative systems with damping matrix C , the dual eigenvalue methodology with linear system solution should be employed to obtain the solution of forced vibration response.
The following two types of systems are discussed: (i) conservative systems without damping, and (ii) nonconservative systems with damping.
i. Conservative systems without damping: Apply the modal analysis for the conservative system to the reduced systems in Eq. (14), without the damping matrix C , to obtain the (n − 1) × (n − 1) modal matrix, U ′ . Proceed with the modal analysis and obtain the forced vibration response (for example, by the inverse Laplace transform, as in the following). Defining q ′ = U ′ η , Eq. (14) becomes The forced vibration response 4 using the inverse Laplace transform is Once the vector of the modal coordinates Finally, the NRB part of the forced vibration response is obtained from the mapping using the constraint matrix in Eq. (11) ii. Nonconservative systems with damping: Following the analytical solution in "NRB motion: free vibration analysis" section, the NRB part of the forced vibration response is the convolution integral in Eq. (31). Note that includes both the generalized coordinates and generalized velocities. Finally, the NRB part of the forced vibration response is obtained from the mapping using the constraint matrix in Eq. (11) Alternative solution using the dual expansion theorem: The dual expansion theorem of the nonconservative eigenvalue systems give rise to Substitute Eq. (33) into the linear system Eq. (29) and pre-multiply by Y T to obtain Apply the orthogonality property to reduce the above equation to Once x(t) is found, the NRB part of the forced vibration response q NRB (t) = S q ′ (t) can be obtained.

Overall response of vibration by combining the RB and NRB solutions
Finally, the solution of the overall response of vibration for such unconstrained system can be determined by combining the appropriate solutions of RB and NRB motions, as follows where q RB (t) represents the solution of the rigid-body motion obtained through formulation presented in "Determine the rigid-body motion (RB)" section, and q NRB (t) is the solution of the non-rigid-body oscillatory motion obtained through formulation presented in "Determine the oscillatory motions (NRB)" section, respectively.

Apply the theory to an example of unconstrained mechanical system
Problem Statement: A three-DoF unconstrained nonconservative mass-spring-damper system is shown in Fig. 2. The generalized coordinates are q = [q 1 q 2 q 3 ] T , where q i are measured from the equilibrium positions. The parameters of the system are: m 1 = 8 kg, m 2 = 2 kg, m 3 = 5 kg , c 2 = 15 N s/m, c 3 = 20 N s/m, c 5 = 10 N s/m , k 2 = 1200 N/m, k 3 = 1500 N/m , and k 5 = 2400 N/m . The initial displacements are q 1 (0) = 0, q 2 (0) = 0.03 and q 3 (0) = 0 m , and the system is released from at rest with zero initial velocity.
Derive the equation of motion for the unconstrained nonconservative system in Fig. 2. Determine (i) the free vibration response and (ii) the forced vibration response subject to a unit impulse Q imp = [δ(t) 0 0] T . Solution: The unconstrained nonconservative system under consideration here is illustrated in Fig. 2.
The equation of motion for the dynamic system shown in Fig. 2 is where q = [q 1 q 2 q 3 ] T is the state vector of the generalized coordinates. It is noted that both the stiffness and damping matrices, K and C , are singular, 5 while the mass matrix, M , is not singular. The system is positive semi-definite. Substituting the given values of the parameters into Eq. (38), we obtain The rigid-body mode is in the null space of the matrix K (to render K u 0 = 0 ); that is, In this case, the rigid-body mode, u 0 , represents a motion in which all three masses move in a synchronized way, as though all three masses were connected as one whole rigid body. The rigid-body mode can be normalized by u T 0 M u 0 . Thus, the normalized rigid-body mode is 6 The corresponding frequency (non-oscillatory) is To solve the eigenvalue problem, we have to reduce the positive semi-definite system by one degree of freedom (from 3 DoF to 2 DoF), by removing the rigid-body mode. Applying the preceding methodology, we can formulate the following constraint equation It is easy to verify that the conventional modal analysis can not be employed in this nonconservative system to obtain the solution of vibration responses. Instead, the general methodology presented in this paper will be employed to obtain the response of free vibration subject to the prescribed initial displacement and velocity.
i. Free vibration response: From Eq. (25), we find Although the system is a nonconservative system, it so happens that u T 0 Cu 0 = 0 . Thus, the free vibration response for the rigid-body motion, β(t) , from Eq. (26) is The response of the rigid-body motion based on Eq. (20) is with the rigid-body (RB) part of the initial conditions being Thus, the initial conditions of the NRB oscillatory motion are Only the first two elements of the NRB initial conditions in Eq. (44) will be used for the initial conditions q ′ (0) and q ′ (0) because q ′ is obtained by eliminating q 3 . 7 The system given in Eqs. (41) and (42) is a 2 × 2 nonconservative positive-definite system, with the initial conditions in Eq. (44). The dual eigenvalue analysis for such nonconservative system will be employed using the methodology in "NRB motion: free vibration analysis" section. Equation (41) is re-arranged according to Eq. (29) to become where x(t) = [q ′ (t)q ′ (t)] T , and the A matrix is 8 The solution of the free vibration response according to Eq. (32) is where the initial condition consists of the first two elements of q ′ and q ′ in q NRB (0) and q NRB (0) because q 3 was eliminated in forming the constraint matrix S.
When q 1 is eliminated instead of q 3 , the initial conditions q ′ (0) and q ′ (0) will take the last two elements of the NRB initial conditions in Eq. (44). When q 2 is eliminated, the initial conditions q ′ (0) and q ′ (0) will take the first and the last elements. 8 The B matrix is zero without excitation in free vibration response.
The solution of free vibration response can be obtained using the x(t) = X e t Y T x(0) , which includes q ′ (t) and q ′ (t) . The response of NRB part of the free vibration can be obtained by q NRB (t) = Sq ′ (t) after simplifying the complex solution, with The overall free vibration response is the sum of the rigid-body (RB) motion in Eq. (43) and the oscillatory motion (NRB) in Eq. (46); that is, q(t) = q RB (t) + q NRB (t) . Therefore, the overall free vibration response is The details of the responses in Eqs. (46) and (47) can be found in Appendix 1. The responses of q(t) in Eq. (47) are linear combination of the RB and NRB responses, and are plotted in Fig. 3. It can be verified that the solution in Eq. (47) satisfies the differential Eq. (38).
Only the general treatment of analysis can be used to obtain the free vibration response. The response is plotted in Fig. 3. The free vibration response plotted in Fig. 3 has a steady-state deviation from zero (or the initial equilibrium position) for the amount of 0.004 m, the same as the rigid-body motion obtained in Eq. (43). The RB motion will cause the entire system to oscillate and asymptotically settles at a new position that is 0.004 m away from the initial equilibrium position. Thus, the free vibration   response of such system resembles a "step response" based on the results shown in Fig. 3.
We have also solved Eq. (38) using the same initial conditions with the numerical differential equation solver NDSolve[] of Mathematica, and overlaid the results by black dashed lines for each DoF in Fig. 3. Both solutions, one from the proposed analytical solution and the other from the numerical solver, are practically the same, which further validates the accuracy of the analytical methodology presented in this paper. The analytical solution, however, provides the physical insights of dynamic characteristics, as revealed in the damping ratios and natural frequencies in the modal space. For example, the low damping ratio, ζ , associated with 1,2 in Table 2 is the main cause of the oscillatory response. ii. Forced vibration response: First, we determine the RB part of the forced vibration response. From Eq. (27), the RB part of the impulse response is  where the impulse excitation is Q(t) = [δ(t) 0 0] T . Note that u T 0 Cu 0 = 0 in this example. Once the response β(t) is found, the RB part of the forced vibration response can be determined by Eq. (20), q RB (t) = u 0 β(t) . Therefore, the RB part of the impulse is Next, let's determine the NRB part of the forced vibration response. Employ the A matrix in Eq. (45) and the right/left eigenvectors in matrices X and Y . The solution of the forced vibration response is given by the convolution integral in Eq. (31). The unit impulse response subject to Q(t) = [δ(t) 0 0] T is The first two elements of x(t) are the reduced state vector q ′ (t) . Thus, the NRB part of the impulse response is q NRB (t) = S q ′ (t) The overall response of forced vibration is the sum of the RB and NRB solutions obtained above. The impulse response is The details of the responses in Eqs. (49), (50) and (51) can be found in Appendix 2.
The unit impulse response is plotted in Fig. 4. Note that the impulse response plotted in Fig. 4 is different from typical impulse responses of a constrained system. Due to the rigid-body motion, the entire system will continue to move away from the initial equilibrium position because of the RB motion in Eq. (48), which is a monotonically increasing motion, while the NRB oscillations are superposed on the RB motion. Thus, the impulse response of a unconstrained system resembles a ramp response of a constrained system.

Application of the theory to impedance control of redundant robots
The theory previously described and exemplified for vibration analysis of multi-DoF mechanical systems can be applied to obtain a closed-form solution of a multi-dimensional joint impedance control system of redundant and non-redundant robotic manipulators. A theoretically sound dynamic response of the robot can be obtained to modulate the behavior of the system through the choice of parameters of the stiffness and damping matrices.

Dynamic response of robotic manipulators performing Cartesian impedance control
In Cartesian impedance control, an n-DoF (redundant or non-redundant) robotic manipulator has a certain level of desired dynamic behavior when interacting with the environment [3]. The equations of motion of an n-DoF redundant manipulator is given by where q is the vector of joint angles, M is the manipulator's mass matrix, G the matrix containing the centrifugal and Coriolis terms and v is the gravity vector term. Choosing τ m as [−K(q)q(t) − C(q)q(t) + v(q) + G(q,q)q(t)] , the system in Eq. (52) becomes: This equation of impedance control is now comparable to Eq. (1) describing a mass, damper and spring system. In order to perform Cartesian impedance control, a desired (52) M(q)q(t) + G(q,q)q(t) + v(q) = τ m + τ ext (53) M(q)q(t) + Cq(t) + Kq(t) = τ ext Fig. 4 Impulse response of the nonconservative system with damping. The response q 1 (t) , q 2 (t) , and q 3 (t) are indicated dynamic behavior is intended in the Cartesian ( m × m ) task space [7] where, K C is the Cartesian stiffness matrix, C C is the damping matrix and M C is the Cartesian mass matrix. The values of K C can be chosen as necessary for the application, and a prospective initial damping matrix C C can be approximated or based on relatively low values. Both Eqs. (53) and (54) assume the similar expression as that in Eq. (1); therefore, the preceding analysis can be employed to obtain dynamic response of impedance control. The impedance control in the joint space is often of primary interests because that is where the physical control is taking place with the motors and encoders. Thus, it is necessary to compute the respective stiffness, damping and mass matrices in the joint space with a specification of Cartesian task requirements.
For a given robot configuration, the inertia mass matrix M in the joint space can be computed [11], and the desired joint stiffness K and damping C matrices can be mapped from the Cartesian task space [5,6,12,13] using the manipulator Jacobian J: where K G = ∂J T ∂q 1 f ∂J T ∂q 2 f · · · ∂J T ∂q n f is the equivalent stiffness due to the change of geometry under the presence of force [5].

Closed-form solution
Similar to the case of mechanical systems, we need to find a solution to Eq. (53) in the joint space. We assume that the external forces f and the changes in the Jacobian J of the manipulator are negligible, and only the first term in Eq. (55) is used to transform the stiffness matrix from Cartesian space to the joint space. With the stiffness and damping, K and C in Eqs. (55) and (56), used in Eq. (1), we can apply Eqs. (2) through (37) in the analysis of joint-space impedance control of robots, with a terminology change outlined in Table 3. The term of "zero-potential-energy, " or the ZP mode, is used in place of the RB motion, with identical physical meaning. The ZP motion will result in no net potential energy change from the stiffness matrix of the impedance control; that is, 1 2 q T ZP K q ZP = 0 , where K is the stiffness matrix of the impedance control in Eq. (53). Likewise, the NZP and NRB motions are also used with identical physical meaning.
This method accounts for the general case of possible redundant and non-redundant manipulators. In the former, this would be equivalent to the unconstrained systems, discussed earlier in the previous section. If the robotic manipulator system is non-redundant, it is no longer unconstrained, the dynamic response would only include the non-zero-potential-energy (NZP) part of the motion. Just like in the case for mechanical systems, the general dynamic response of the robot can be expressed as a linear combination of all the bases, u i , as in Eq. (6).

ZP and NZP motions
The Given that in general both K and C matrices go through similar transformations via the Jacobian matrix J, the second term in Eq. (58) will vanish. When u 0 is normalized, this equation can be simplified as β (t) = u T 0 τ ext , and solved by including the initial conditions in this reduced space. To determine the oscillatory motions, we remove the known ZP mode q i by imposing upon the state vector q(t) a n×(n − 1) constraint matrix, S , as introduced earlier from Eq. (6) to Eq. (19), to build the M ′ , C ′ and K ′ matrices.
Once the NZP system is positive definite, we can obtain the NZP responses as described in the theory for mechanical discrete systems. After that, it is necessary to combine both ZP and NZP solutions, similar to Eq. (6) to finally obtain the full dynamic response. Then, a (57) q ZP (t) = u 0 β(t) proper Cartesian damping matrix can be chosen [4]. Note that the complex conjugate eigenvalues render the natural frequency and damping ratio in the modal space, as follows where a and b are the real and imaginary parts of the eigenvalues obtained in Eq. (31), and γ = tan −1 (b/a) . The dynamic response of the robotic system performing Cartesian or joint-based impedance control is determined by the damping ratios, ζ , and the damped frequency, ω d , through a and b.

Special case of multiple redundancies
For cases when r = n − m > 1 , the methodology previously presented can also be applied; however, it is necessary to successively remove all r ZP modes from the system, while satisfying the orthogonality property with the mass matrix M , as in Eq. (7), where each u i represents a ZP mode.
With multiple redundancy and constraint matrices from the ZP modes, we have The r ZP modes of the system become an (n × r) matrix U r containing the vectors of the null space of the stiffness matrix K . There are two main ways of obtaining the NZP motions, as follows.

Successive removal of every ZP mode
Each ZP mode is removed one at a time, reducing the redundancy by one each time. This requires building a set of matrices using the constraint matrix each time until all r ZP modes are removed and the system becomes positive definite, as in the following.

Abbreviated removal of ZP modes
Alternatively, all ZP modes can be removed by using the full U r matrix to reduce the system directly into a (n − r) × (n − r) positive definite system. The following equations explain the procedures.
where u 1 to u r−1 are derived with the following state variable transformations Both ways of removing ZP modes are further explained in the following example of robotic impedance control.

Example of a redundant 7-DoF robot
The example used to illustrate the analysis is a general robotic system which consists of 7 DoF. We choose to consider a Cartesian impedance control with 5 DoF to demonstrate the procedures to handle multiple redundancy. Therefore, we have n = 7 and m = 5 in this case, giving two-DoF redundancy with r = 2 . Both methods of successive removal of the ZP modes in (i) and the abbreviated method in (ii) are employed.
At a chosen robotic configuration, there is a Jacobian matrix J associated with it, as well as a mass matrix M , and a chosen set of stiffness and damping matrices K and C mapped from the Cartesian space. All these matrices can be found in Appendix 3, at the end of the paper. The 7 × 7 impedance control equation in the joint space is where the mass, damping and stiffness matrices are:

Mq(t) + Cq(t) + K q(t) = Q
As stated in the theory, the characteristic of an unconstrained, positive semi-definite system is a singular stiffness matrix, K . For this case, the eigenvalues of the stiffness matrix are: This system is not solvable, unless the two ZP modes (r = 2) are removed. Once this is done, we can solve for the oscillatory and non-oscillatory motions as presented in Eqs.

Successive removal of every ZP mode
The stiffness matrix K is associated with two ZP modes. Both ZP modes belong to the null space of K . When normalized and making them orthogonal to each other, we obtain In order to remove the first ZP mode we use the first vector u 0 to build its corresponding constraint matrix S 0 , such that With this constraint matrix S 0 we can now build a new system free of this first ZP mode: The numerical values for these matrices are in Appendix 3. This 6 × 6 system in the reduced space q ′ is still not free from the second ZP mode u 1 . The system is not positive definite, with the eigenvalues of In order to remove the remaining ZP mode, we need to build a second constraint matrix S 1 by first obtaining and normalizing this ZP mode from K ′ in the next reduced space of K ′ u ′ 1 = 0 and S 0 u ′ 1 = u 1 . Once this ZP mode is obtained we can build S 1 with With this constraint matrix S 1 the new system in the q ′′ space can be found The numerical values of these matrices are in Appendix 3. Now, this system is positive definite since both ZP modes have been removed, with the eigenvalues of This ( 5×5 ) system is now solvable by applying the theory of mechanical vibration of discrete systems presented ealier.
Next, we are going to show that the abbreviated method is equivalent to removing the ZP modes successively one at a time.

Alternative method to remove ZP modes
From Eqs. (61) and (62) and information of every ZP mode, we can build the constraint matrices. For S 0 , we can directly use u 0 We can now use S 0 to obtain the next ZP mode in the reduced space q ′ with u ′ 1 = S * 0 u 1 , where S * 0 is the pseudo-inverse of the first constraint matrix, and u 1 is the next ZP mode in the original ( 7 × 7 ) space. It is important to always work with the ZP vectors normalized with respect to M , and the inverse of S 0 must be unique and consistent with M , i.e. S * 0 = (S T 0 MS 0 ) −1 S T 0 M . Once this ( 6 × 1 ) ZP mode, u ′ 1 , is obtained, the S 1 constraint matrix, with q ′ = S 1 q ′′ , can be constructed as follows This is exactly the same as the S 1 obtained in Eq. (63) from successive removal of the ZP modes. The full constraint matrix S can be obtained without needing to build systems in each step, and as in Eq. (60) With S = S 0 S 1 , the matrix can be written as With this full constraint matrix S , we can obtain the positive definite system directly as follows The results are identical to the system obtained by using the method of building a new system by successively removing each ZP presented earlier. This reduced system leads to the following damping ratios and natural frequencies. It is noted that there are 3 modes with small damping ratios, ζ , in Table 4.

Solve the positive definite system
This positive definite system {M ′′ , C ′′ , K ′′ } now can be solved by using Eqs. (29) to (32) of the dual eigenvalue analysis. An initial displacement of the following is given to the system The response of NZP motion q ′′ 1 to q ′′ 5 denoted by can be obtained by following the methodology of the dual eigenvector analysis. The results of q ′′ 1 to q ′′ 5 are in Appendix 3.
The dynamic response of the physical system q(t) can now be obtained by the following equation where the NZP component from Eq. (37) is q NZP = S 0 S 1 q ′′ (t) . The results are plotted in Fig. 5 where both ZP and NZP parts of the response are combined in the responses from q 1 to q 7 .
The dynamic response in Fig. 5 is plotted for the first five seconds to illustrate the characteristics of free vibration response. The results are indicative of the typical  behavior and can give a glimpse of the characteristics of dynamic response of the impedance control through the overshoots and settling times of the response. When plotted over a longer period of time, we can examine the low damping characteristics with low frequency that takes longer to settle.

Discussion
Free vibration responses of a unconstrained and nonconservative discrete mechanical system with prescribed initial displacement show that the responses resemble a step response of a system with external force, due primarily to the RB mode of the unconstrained system. As shown in the plots of the combined RB and NRB motions, the initial conditions are satisfied. Likewise, an impulse response of a unconstrained mechanical system has a forced vibration response that resembles a ramp response of forced vibration, due to the RB mode. It is noted that the coefficient of the β term in Eq. (23) will not be unity unless the RB mode, u 0 , is normalized with respect to the mass matrix, as defined in Eq. (7). In the case of free vibration, when u 0 is not normalized according to Eq. (7), the non-unity coefficient will disappear from the homogenous equation and can result in incorrect RB component of the initial condition. On the other hand, forced vibration in Eq. (22) will not be affected because the coefficient will be carried through with the non-zero forcing function term u T 0 Q. In applying the methodology to an example of impedance control of a Baxter robot, the complex eigenvalues of the 5×5 positive definite system in the reduced space of robotic impedance control are −1.450 ± i14.54 , −0.3116 ± i3.700 , −0.1833 ± i3.066 , −0.3660 ± i2.046 and −0.04169 ± i0.8158 . These complex eigenvalues result in damping ratios, ζ , and natural frequencies, ω n listed in Table 4. It is noted that the damping ratios are low, resulting in highly oscillatory responses, as shown in the responses plotted in Fig. 5. Since the responses of q 1 to q 7 are the linear combination of the bases obtained from the 5×5 positive definite system, the lower damping ratios, especially when accompanied by low natural frequency, will appear in q 1 to q 7 to dominate the response. It is important to be able to increase the damping ratios in order to improve the overall dynamic responses.
This can be achieved, with a given stiffness matrix, by choosing different values in the elements of damping matrix. However, study shows that increasing all elements of the damping matrix may not necessarily improve all damping ratios. Furthermore, certain elements of the damping matrix have more influence on certain damping ratios of the complex eigenvalues. The methodology presented in this paper, nevertheless, can assist in determining which diagonal and/or off-diagonal elements in the damping matrix can improve the damping ratios. This reduces blind guessing or trial-and-error in choosing values of the elements of the damping matrix [4,6] to improve dynamic response.
By applying the dual eigenvalue analysis in Eqs. (29) to (32), we are able to evaluate the dynamic response of the impedance control by criteria such as damping ratios and frequency of the NZP part of the motion. This methodology with the closed-form solution allows us to adjust the stiffness and damping parameters in the K and C matrices, respectively, according to the required applications of the performance of dynamic responses. It also allows us to modulate the dynamic response of the null space of the robotic task [14][15][16]. If there are any desired additional tasks for the redundant DoF(s), the stiffness and damping parameters can be mapped into the joint space by the use of a null space projection matrix, and from there on, proceed as shown here. In summary, the methodology presented in this paper can provide the advantage to systematically adjust the elements in K and C matrices of impedance control to modulate the dynamic responses of the system.
As a case in point, when trying to perform impedance control in the Cartesian ( 6 × 6 ) space using a redundant 7-DoF manipulator at a chosen configuration, the given stiffness for the task is An initial damping matrix is picked as The following damping ratios can be obtained: 0.0108, 0. 0039, 0.0162, 0.0342, 0.116 and 0.0953. In contrast, if we apply the methodology presented here to select a damping matrix, one possible non-diagnoal damping matrix, after a synthesis using the methodology, is The damping ratios of C C1 are: 0.0326, 0.1081, 0.1335, 0. 3419, 0.1342 and 0.1274. This revised impedance control with the new C C1 damping matrix has a better dynamic response because of the improvement of damping ratios of the multiple DoF system.

Conclusion
A general methodology to obtain a closed-form solution for free and forced vibration responses are presented for unconstrained mechanical system and robotic manipulators under impedance control. The methodology includes a novel procedure to deal with nonconservative unconstrained systems or redundant robots containing rigid body (RB) or zero potential energy (ZP) modes, which must be removed, using constraint matrix, to obtain the dynamic response. After that, the non-positive definite system can be mapped via constraint matrix to form a positive definite system in a reduced space, to obtain the complete dynamic response of the system.
When applying to robotic manipulators performing Cartesian and joint-based impedance control, this methodology can provide the analysis and synthesis of dynamic characteristics of the impedance control system. This theory can help in choosing the elements of the damping matrix to attain prescribed dynamic response with damping characteristics without guessing or trialand-error. The elimination of trial-and-error in determining the elements of the damping matrix in impedance control to meet the requirements of dynamic response, by using this novel analytical methodology, is another important contribution of this paper. The methodology was illustrated by examples of both a discrete mechanical system and a redundant robotic manipulator in impedance control with damping modulation.