## 2 Description of Dynamic Model for a Robot with Elastic Joints  
Under the assumptions in [PERSON], [PERSON] et al. [29], the dynamic model of a robot with elastic joints can be expressed as  
\[\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}}+\mathbf{C}\mathbf{\left(\mathbf{q},\dot{\mathbf{q}}\right)}+\mathbf{G} \mathbf{\left(\mathbf{q}\right)}+\mathbf{K}_{\mathbf{q}}\mathbf{\left(\mathbf{q}-\mathbf{\theta}\right)}+ \mathbf{F}_{\mathbf{q}}\dot{\mathbf{q}}+\mathbf{D}\mathbf{\left(\dot{\mathbf{q}}-\dot{\mathbf{\theta}} \right)}=-\mathbf{J}^{T}\mathbf{\left(\mathbf{q}\right)}\mathbf{f}, \tag{1}\]  
\[\mathbf{B}\ddot{\mathbf{\theta}}+\mathbf{K}_{\mathbf{q}}\mathbf{\left(\mathbf{\theta}-\mathbf{q}\right)}+ \mathbf{F}_{\mathbf{\theta}}\dot{\mathbf{\theta}}+\mathbf{D}\mathbf{\left(\dot{\mathbf{\theta}}-\dot{ \mathbf{q}}\right)}=\mathbf{\tau}, \tag{2}\]  
where \(\mathbf{q}\), \(\dot{\mathbf{q}}\), and \(\ddot{\mathbf{q}}\) are the \(\left(n\times 1\right)\) vectors of position, velocity, and acceleration of links, respectively. \(\mathbf{\theta}\), \(\dot{\mathbf{\theta}}\), and \(\dot{\mathbf{\theta}}\) are the \(\left(n\times 1\right)\) vectors of positions, velocities, and accelerations of motors. \(\mathbf{M}(\mathbf{q})\) is the \(\left(n\times n\right)\) robot link inertia matrix, \(\mathbf{C}\mathbf{\left(\dot{\mathbf{q}},\dot{\mathbf{q}}\right)}\) is the \(\left(n\times 1\right)\) vector of centrifugal and Coriolis torques, \(\mathbf{K}_{\mathbf{q}}\) is the \(\left(n\times n\right)\) diagonal matrix of joint stiffness coefficients, \(\mathbf{G}\mathbf{\left(\mathbf{q}\right)}\) is \(\left(n\times 1\right)\) vector of gravitational torque, and \(\mathbf{B}\) is a \(\left(n\times n\right)\) constant diagonal matrix, including the rotor inertia in the gear ratios. \(\mathbf{F}_{\mathbf{q}}\), \(\mathbf{F}_{\mathbf{\theta}}\), and \(\mathbf{D}\) contain, in this order, the viscous coefficient in the link side, the viscous coefficient in the motor side, and the damping of the elastic springs at the joints. \(\mathbf{\tau}\) is the \(\left(n\times 1\right)\) input vector of driving torques, \(\mathbf{f}\) is the \(\left(n\times 1\right)\) vector of contact force exerted by the end effector on the environment, and \(\mathbf{J}\mathbf{\left(\mathbf{q}\right)}\) is the Jacobian matrix that relates joint velocities \(\dot{\mathbf{q}}\) with the vector of end-effector velocities, \(\dot{\mathbf{x}}\). The Jacobian is assumed to be non-singular. The transpose of the matrix, \(\mathbf{J}^{T}\mathbf{\left(\mathbf{q}\right)}\), also relates the end effector force with the joint torques.  
\[\dot{\mathbf{x}}=\mathbf{J}\mathbf{\left(\mathbf{q}\right)}\dot{\mathbf{q}}, \tag{3}\]  
\[\mathbf{\tau}=\mathbf{J}^{T}\mathbf{\left(\mathbf{q}\right)}\mathbf{f}, \tag{4}\]  
For analysis purposes, the environment is modelled as a frictionless and elastically compliant plane, which is very common in force control [18; 29; 30]. One contact point is considered, and the contact force is expressed as  
\[\mathbf{f}=\mathbf{K}_{\mathbf{e}}\mathbf{\left(\mathbf{x}-\mathbf{x}_{\mathbf{e}}\right)}, \tag{5}\]  
where \(\mathbf{x}\) is the end-effector position, \(\mathbf{x}_{\mathbf{e}}\) is the position at the contact point, and \(\mathbf{K}_{\mathbf{e}}\) is the constant symmetric stiffness matrix of the environment.