Student Member, IEEE Marco Faroni

Manuel Beschi

Nicola Pedrocchi

Antonio Visioli

Predictive Inverse Kinematics for Redundant Manipulators with Task Scaling and Kinematic Constraints

Keywords: Inverse kinematics, trajectory scaling, redundant robots, task priority, model predictive control, motion planning

I. INTRODUCTION

The inverse kinematics (IK) of robot manipulators permits to determine the joint variables associated with a given task (e.g., in the Cartesian space). Generally, the solution of this problem is not trivial, since it involves the inversion of multivariable nonlinear functions. Moreover, when the robot is kinematically redundant, the problem becomes even more complex: the system is underdetermined and it admits infinite solutions. Among them, the IK algorithm should choose the one that tends to satisfy secondary tasks or, similarly, optimize a given performance index (e.g., the manipulability). Typically, the IK is solved at velocity or acceleration level. As a matter of fact, at differential level, it is possible to exploit the linear transformation between joint and task velocities/accelerations given by the Jacobian. Then, the corresponding joint configuration is calculated by forward integration. The same occurs for redundant manipulators. In this case, the joint variables are obtained by means of a generalized inverse (usually the Moore-Penrose pseudoinverse) of the Jacobian, which is in this case a rectangular matrix with more columns than rows. A secondary task is usually handled by the projection in the null space of the Jacobian, in order not to interfere with the execution of the primary task [START_REF] Nakamura | Advanced Robotics: Redundancy and Optimization[END_REF], [START_REF] Liegeois | Automatic supervisory control of the configuration and behavior of multibody mechanisms[END_REF]. This formulation has been generalized for any number of tasks in [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF], giving birth to the so-called task-priority framework. A great disadvantage of this method is that it does not take explicitly into account the physical limits of the robot. Thus, the obtained solution may violate such limits, resulting infeasible for the real capabilities of the robot. In the case of redundant manipulators, this issue has often been handled by converting the hard bounds into soft ones in a cost function that tends to keep the joints in the middle of their ranges [START_REF] Liegeois | Automatic supervisory control of the configuration and behavior of multibody mechanisms[END_REF], [START_REF] Atawnih | Kinematic control of redundant robots with guaranteed joint limit avoidance[END_REF] or by using a weighted pseudo-inverse to limit velocities or accelerations [START_REF] Chan | A weighted least-norm solution based scheme for avoiding joint limits for redundant joint manipulators[END_REF]. However, these methods cannot ensure the satisfaction of the bounds, and therefore they do not prevent from having infeasible solutions. Note that, in the practical case, this would lead to the saturation of the actuators, with consequent deformation of the original task or even the stop of its execution.

Recent works cope with this by addressing the inverse kinematics as a constrained optimization problem [START_REF] Rocchi | Opensot: a whole-body control library for the compliant humanoid robot Coman[END_REF]. In general, the task is converted into a Least Square Problem (LSP) and the solution is searched inside an admissible set. Secondary tasks are handled similarly, by successive LSP, wherein the admissible set is limited by the solutions of the tasks with higher priority. A general formulation M. Faroni and A. Visioli are with the Dipartimento di Ingegneria Meccanica e Industriale, Università degli Studi di Brescia, Brescia, 25123 Brescia, Italy (e-mail: m.faroni003@unibs.it; antonio.visioli@unibs.it).

M. Beschi and N. Pedrocchi are with Istituto di Sistemi e Tecnologie Industriali Intelligenti per il Manifatturiero Avanzato, Consiglio Nazionale delle Ricerche, 20133 Milan, Italy (e-mail: manuel.beschi@stiima.cnr.it; nicola.pedrocchi@stiima.cnr.it).

following this paradigm is presented in [START_REF] Kanoun | Kinematic control of redundant manipulators: generalizing the task priority framework to inequality tasks[END_REF], where the priority framework is extended to inequality tasks as well. A great advantage of the LSP formulation is that it can explicitly include constraints in a Quadratic Program (QP), which can be solved efficiently by numerical QP solvers (e.g., [START_REF] Ferreau | qpOASES: A parametric active-set algorithm for quadratic programming[END_REF]- [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF]). In this way, the solution (if it exists) will surely satisfy the hard bounds of the manipulator. However, the use of such approach requires particular attention on the way the constraints are posed. In particular, if joint configuration, velocity and acceleration limits are considered, viability of the states must be ensured. The viability of a set of states implies that, given a state within the viable set, all the future states will belong to that set [START_REF] Blanchini | Set invariance in control[END_REF]. This is fundamental to ensure that constraints will be satisfied in the future. This problem and issues related to its discrete-time implementation have been discussed in details in [START_REF] Del Prete | Joint position and velocity bounds in discrete-time acceleration/torque control of robot manipulators[END_REF].

The level of satisfaction of the tasks (included the primary one) also strongly depends on how the task is defined with respect to the constraints. In other words, a too-demanding task would rapidly lead to the saturation of the joint commands, with consequent considerable deformation of the original task. Note that the primary task imposed to a manipulator typically consists in following a predefined path in the Cartesian space. Therefore, a significant deformation of such task is often unacceptable and should be avoided by the IK algorithm by any means. A widespread strategy to partially overcome this issue consists in the time scaling of the trajectory. This allows the time dilation of the task, in order to fulfill both the geometrical task and the joint bounds. Several task-scaling techniques have been proposed in the literature. First examples appeared in [START_REF] Hollerbach | Dynamic scaling of manipulator trajectories[END_REF], [START_REF] Dahl | Torque-limited path following by online trajectory time scaling[END_REF], in order to cope with torque boundaries of the actuators in a dynamic-based control of a manipulator. Another approach was proposed in [START_REF] Gerelli | Nonlinear variable structure filter for the online trajectory scaling[END_REF], [START_REF] Lo Bianco | Online trajectory scaling for manipulators subject to high-order kinematic and dynamic constraints[END_REF], where the problem was reformulated in terms of nonlinear filtering of the joint trajectories, based on kinematic and/or dynamic constraints. As regards pure kinematic control, an important work is represented by [START_REF] Antonelli | A new on-line algorithm for inverse kinematics of robot manipulators ensuring path tracking capability under joint limits[END_REF], where a scaling factor is introduced in the computation of the IK so that the joint commands can be scaled within the boundaries. Nonetheless, if the robot is redundant, this solution might be nonoptimal, as other velocities/accelerations (among the infinite possible solutions) could allow a better or even full satisfaction of the task. An improved method was proposed in [START_REF] Flacco | Control of redundant robots under hard joint constraints: Saturation in the null space[END_REF], wherein an efficient exploration of the null space permits to find the solution ensuring the best scaling factor. This problem is equivalent to a QP where the task is expressed as an equality constraint, which forces the executed task to follow the same direction of the desired one, whereas the scaling factor is optimized in the cost function. The same algorithm was developed at acceleration level in [START_REF]Motion control of redundant robots under joints constraints: Saturation in the null space[END_REF], in order to recover continuity on the velocity and to include acceleration limits. However, the introduction of acceleration bounds does not always guarantee the existence of a feasible solution. In this case, the equality constraint has to be converted again into an LSP, allowing the modification of the primary task. Therefore, the obtained solution will firstly minimize the deformation with respect to the original geometrical path, and secondly will optimize the scaling factor.

The methods illustrated up to this point are all termed as local, since they only take into account information available at the current iteration. This means that the solution is optimal with respect to the current configuration of the manipulator, but not with respect to the global execution of the task. As a matter of fact, the iterative local resolution of the inverse kinematics might move the robot to a disadvantageous configuration with respect to the global task, as shown in [START_REF] Faroni | A predictive approach to redundancy resolution for robot manipulators[END_REF]. Similarly, the local optimization of the task scaling might lead to worse results if compared to an optimization involving the future evolution of the task, the constraints and the kinematic variables. Anyhow, the use of global methods for the redundancy resolution or the scaling problem (i.e., the modification of the timing law) is strongly limited by the heavy computational load, which makes them unsuitable for online implementation.

To combine the strengths of local and global algorithms, this paper proposes a unified method to solve the inverse kinematics and the task scaling problem of a redundant manipulator by means of a model predictive approach. Model Predictive Control (MPC) is becoming more and more popular in robotic manipulation and mechatronics, thanks to the growing computational power of modern processors [START_REF] Wilson | Nonlinear model predictive control schemes with application on a 2 link vertical robot manipulator[END_REF]. In this context, contributions usually address low-level control problems of actuators for trajectory tracking [START_REF] Hedjar | Nonlinear Receding-Horizon control of rigid link robot manipulators[END_REF]- [START_REF] Ferrara | A robust MPC/ISM hierarchical multi-loop control scheme for robot manipulators[END_REF]. Few works also apply predictive techniques to high-level problems such as pointto-point trajectory generation [START_REF] Ghazaei Ardakani | Real-time trajectory generation using model predictive control[END_REF] or motion planning for mobile manipulators [START_REF] Buizza Avanzini | Reactive constrained model predictive control for redundant mobile manipulators[END_REF].

In general, MPC is based on the optimization of a cost function over a finite predictive horizon, in order to get an optimal set of control actions for the considered time interval. Then, the first control action is applied to the system and the procedure is repeated iteratively, in the spirit of receding horizon control. One great advantage of MPC techniques is the efficient handling of constraints. Moreover, the predictive approach permits to take into account the evolution of the system in the successive time instants as well. For this reason, such approach seems to be a good candidate for overcoming the problems of local scaling and punctual resolution of the IK in a redundant manipulator. Following this principle, [START_REF] Schuetz | Predictive online inverse kinematics for redundant manipulators[END_REF] exploited Pontryagin's Minimum Principle to set up a nonlinear optimization problem over a finite horizon. Nonetheless, kinematic bounds are not taken into account and this represents a considerable limitation of the algorithm. As a matter of fact, their inclusion would lead to a much heavier optimization problem, increasing significantly the computational time. In [START_REF] Tassa | Control-limited differential dynamic programming[END_REF], [START_REF] Koenemann | Whole-body model-predictive control applied to the HRP-2 humanoid[END_REF], nonlinear MPC was applied to the operational-space inverse dynamics problem of a humanoid robot with bounded torques. The algorithm exploits differential dynamic programming to set up a nonlinear problem and the use of a fast nonlinear solver and powerful hardware permit to reach a computational time of 50 ms. In [START_REF] Faroni | A predictive approach to redundancy resolution for robot manipulators[END_REF], a fast predictive approach to redundancy resolution of robot manipulators was proposed: the task was considered as a linearized constraint and an ad-hoc continuoustime formulation of the predictive equations was derived to lighten the computational burden of the resulting QP. Although the final QP is an approximation of the original nonlinear problem, it allows the computational time to be of the same order of magnitude of local methods (in the order of 1 ms), but still conferring the advantages of a predictive strategy. However, when the robot is required to perform demanding trajectories that often lead to joint saturations, this strategy might provoke a big deformation of the original task. In particular, if infeasible motions are required to the robot, this approach is not able to adapt the desired trajectory to the real capabilities of the robot.

For this reason, in this paper we propose a predictive scaling technique to preserve the geometrical task when the desired motion is not realizable by the robot. In this sense, the proposed method confers to the IK technique the advantages of both a predictive redundancy resolution and a predictive scaling. Then, a Closed Loop Inverse Kinematics (CLIK) scheme is implemented [START_REF] Chiacchio | Closedloop inverse kinematics schemes for constrained redundant manipulators with task space augmentation and task priority strategy[END_REF], in order to recover from position errors (given by numerical integration or out-of-path initial configurations). Moreover, a viability analysis is conducted and included in the online QP as a linear constraint. In this way, the robot states are always ensured to stay within a feasible region. This, along with considering the task in the cost function (and not as a hard constraint), always allows the feasibility of the QP.

Compared to [START_REF] Faroni | A predictive approach to redundancy resolution for robot manipulators[END_REF], the proposed technique permits to obtain smaller deformation of the geometrical path, especially when the desired task is demanding with respect to the robot capabilities. Furthermore, compared to local scaling algorithms, it is possible to obtain smaller deformations of the desired tasks and higher mean values of the scaling factor as well.

II. PRELIMINARIES

A. Inverse kinematics as an optimization problem

Consider an n-degree-of-freedom manipulator and an mdimensional task x with m ≤ n, defined with respect to a curvilinear abscissa σ , corresponding to the nominal trajectory time. In order to perform the specified task, the following equation must hold for all values of σ :

x(σ

) = f (q(σ )) (1)
where x ∈ R m is the task position vector, q ∈ R n is the joint configuration vector, and f is the forward kinematics of the robot.

Differentiating with respect to σ gives:

x (σ ) = J(q) q (σ ) (2)
where x = dx/dσ , q = dq/dσ , and J = ∂ f /∂ q is the m × n Jacobian matrix. Since σ corresponds to the nominal desired time, then x represents the desired task velocity, and q corresponds to the joint velocity. The IK problem consists in finding the joint velocities q that satisfies the linear system (2). This can be reformulated as an LSP problem such that:

q = argmin q ∈S 1 2 J q -x 2 ( 3
)
where S is an admissible set of joint velocities, which can include the kinematic boundaries of the manipulator. Robot limits at the current iteration are typically linearized around the current state of the system, as described in [START_REF] Cheng | Resolving manipulator redundancy under inequality constraints[END_REF]. In this way, S results in a set of linear constraints and (3) therefore results to be a QP. When the robot is redundant with respect to the first task, a secondary task E q = g can be specified. The new optimization problem is:

q = argmin q ∈S 2 1 2 E q -g 2 ; S 2 = argmin q ∈S 1 2 J q -x 2 . ( 4
)
This approach can be repeated iteratively for a general number of tasks [START_REF] Kanoun | Kinematic control of redundant manipulators: generalizing the task priority framework to inequality tasks[END_REF]. The prime goal of an efficient IK algorithm is the maximal satisfaction of the primary task, and this depends on the solution of (3). In particular, if the kinematic limits specified in S lead to a constrained solution of the LSP, a spatial deformation of the task arises. To reduce this phenomenon, a modification of the timing law might be performed. As proposed in [START_REF] Flacco | Control of redundant robots under hard joint constraints: Saturation in the null space[END_REF], at each time instant, (2) can be modified as:

s x (σ ) = J(q) q(t) (5)
where s = dσ /dt can be seen as a scaling factor that deforms the velocity profile but without modifying the path. Since only downscaling is considered, the value of s has to be between 0 and 1. When no scaling occurs, s is equal to 1, whereas s smaller than 1 slows down the original trajectory. The desired task velocity x is calculated over the parameter σ = s dt. Since in the discrete-time implementation of the algorithm, s is calculated as a scalar value at each time instant, the parameter σ is computed at each cycle as

σ (k) = σ (k -1) + T s(k -1)
, with σ (0) = 0, s(0) = 1 and T is the sampling period of the controller. In this way, the primary objective of the IK problem becomes the satisfaction of the geometrical path of the desired task. Secondly, the scaling factor will be optimized, in order to maximize the fulfillment of the original primary task. Finally, the secondary task will be accomplished as well as possible. This hierarchy of QPs can be formalized and extended to a generic number of secondary tasks in terms of a lexicographic multi-objective quadratic problem [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF]. In lexicographic optimization, an ordered sequence of cost functions is optimized with priority determined by the order of the sequence. By defining an ordered sequence {w 1 , w 2 , w 3 , . . . , w h }, h ∈ N, where:

w 1 := 1 2 J q -sx 2 , (6)
w 2 := 1 2 (1 -s) 2 , (7)
w i := 1 2 E i q -g i 2 i = 3 . . . h, (8)
the terms w i , i = 3, . . . , h represent generic secondary tasks with decreasing priority. The scaling problem can be formalized as follows:

lex min q, s {w 1 ,w 2 , w 3 , . . . , w h } subject to q ∈ S , s ∈ [0, 1]. ( 9
)
This approach is equivalent to the one developed in [START_REF] Flacco | Control of redundant robots under hard joint constraints: Saturation in the null space[END_REF]. In that case, the higher-priority task could be substituted by the equality constraint J q = sx , as acceleration limits were not considered.

In [START_REF]Motion control of redundant robots under joints constraints: Saturation in the null space[END_REF] the algorithm was formulated at acceleration level as well, in order to recover continuity on the joint velocities and to include acceleration bounds. However, in this way, feasibility is not ensured anymore, as explained in the following section.

B. Definition of the constraints and existence of the solution

he kinematic variables of the robot are typically required to respect the joint configuration and velocity limits, that is:

Q min ≤ q ≤ Q max (10) Qmin ≤ q ≤ Qmax ( 11
)
where Q min , Q max ∈ R n are the minimum and maximum joint configurations respectively, while Qmin , Qmax ∈ R n are the minimum and maximum joint velocities respectively. Moreover, if acceleration limits are taken into account, it must be:

Qmin ≤ q ≤ Qmax ( 12
)
where Qmin , Qmax ∈ R n are the minimum and maximum joint accelerations respectively. However, to ensure that each joint can stop within its configuration bound without exceeding the maximum acceleration Qmax,i , the following viability condition must be imposed (see [START_REF] Scheint | Invariance control in robotic applications: Trajectory supervision and haptic rendering[END_REF], [START_REF] Zanchettin | Motion planning for robotic manipulators using robust constrained control[END_REF]):

Ψ(q i , qi ) ≤ 0 (13)
where, for i = 1 . . . n (i.e., for all joints),

Ψ(q i , qi ) :=    q i + q2 i 2 Qmax,i -Q max,i if qi > 0 q i -Q max,i otherwise
and likewise for the lower bound. Issues in the discrete-time implementation of such constraints are well detailed in [START_REF] Del Prete | Joint position and velocity bounds in discrete-time acceleration/torque control of robot manipulators[END_REF]. Note that the use of a scaling factor as in [START_REF] Flacco | Control of redundant robots under hard joint constraints: Saturation in the null space[END_REF] always guarantees the existence of a solution to the task J q = sx (at worst q = 0 and s = 0), since acceleration bounds were not considered. Therefore, assuming J is full rank, the task could be implemented as an equality constraint in the QP problem, ensuring the full satisfaction of the geometrical task. Nonetheless, there is no control on the acceleration values and this means the algorithm even admits discontinuities between successive solutions. When also acceleration limits are taken into account, the presence of a scaling factor does not ensure the existence of a solution that satisfies both the kinematics limits and the geometrical task. Therefore, a LSP solution would be needed again, with consequent possible deformation of the desired path.

Note that bounded accelerations play an important role in the suitability of planning algorithms. As a matter of facts, acceleration limits are often exploited in kinematic control to approximate torque limits (which typically represent the physical dynamic limits of the actuators). Generally, approximated (and conservative) constant limits for the acceleration are imposed, based on some estimation rules (e.g., the one proposed in [START_REF] Del Prete | Joint position and velocity bounds in discrete-time acceleration/torque control of robot manipulators[END_REF]), as in [START_REF] Antonelli | A new on-line algorithm for inverse kinematics of robot manipulators ensuring path tracking capability under joint limits[END_REF], [START_REF]Motion control of redundant robots under joints constraints: Saturation in the null space[END_REF], [START_REF] Scheint | Invariance control in robotic applications: Trajectory supervision and haptic rendering[END_REF].

III. METHOD

Purely local methods are based on the evaluation of the current configuration of the system. In this section, we reformulate the IK problem in a model-predictive-like framework, in order to take into account also the evolution of the variables along the predictive horizon. The basic idea of the method consists in formulating the resolution of the predictive IK problem as a constrained QP, which has to be solved at each time step. Then, in the spirit of receding horizon control, the first input of the sequence is applied to the system and the procedure is repeated until the end of the trajectory.

A. Formulation of the predictive equations

A linear model of the joint-space kinematic variables for an an n-degree-of-freedom manipulator can be derived by considering each joint as a r-th order integrator and, then, by coupling all the joints in state-space representation. Then, a linear predictive model can be easily derived by forward solving the resulting difference equation, as in classic linear MPC. Due to small sampling periods typical of robotic systems, sufficiently long horizons with uniform sampling periods would lead to a very large number of control variables involved in the predictive model [START_REF] Dimitrov | On the implementation of model predictive control for on-line walking pattern generation[END_REF]. In order to reduce the complexity of the predictive model, the complexity reduction strategy proposed in [START_REF] Faroni | A predictive approach to redundancy resolution for robot manipulators[END_REF] has been applied. Such technique uses just a small number of prediction and control time instants along the horizons, instead of using all the sampling times. In particular, given an ordered set of prediction time instants {τ i }, i = 1, . . . , p, and an ordered set of control time instants {t j }, j = 1, . . . , c, the predictive horizon (in terms of seconds) corresponds to τ p , and the control horizon (in terms of seconds) corresponds to tc . The order r of the integrators determines whether the control input is velocity (r = 1), acceleration (r = 2), and so on. Then, the predictive equations at the current time instant t 0 can be written as:

q (d) (τ) = L d z 0 + F d u (t) d = 0, . . . , r -1 (14)
where q

(d) (τ) ∈ R np×1
is the prediction of the d-th derivatives of the joint configuration at times t 0 + τ i , z 0 is the state vector at time t 0 , being

z := [q 1 , q1 , . . . q (r-1) 1 , . . . , q n , qn , . . . q (r-1) n ] T ∈ R n(r-1)×1 , ( 15
)
u (t) ∈ R nc×1
is the input function (i.e., the r-th derivative of the joint configuration) at times t 0 + t j , and L d ∈ R np×nr and F d ∈ R np×nc are predictive matrices derived as in [START_REF] Faroni | A predictive approach to redundancy resolution for robot manipulators[END_REF].

B. Choice of the predictive and control time instants

The predictive model developed in the previous section is based on the a-priori placement of the prediction and control time-instants τ i and t j along the prediction and the control horizons respectively. Although prediction and control time-instants can be chosen to be different, hereafter they have been considered to be equal without lose of generality. Since the reference signal (i.e., the Cartesian task) is, in general, non-repetitive, it is not possible to determine an optimal placement of τ i and t j . For this reason we use a simple empiric tuning rule based on the following guidelines:

• the control and prediction time-instants have the same distribution (i.e., τ i = ti ); • the first time-instant is equal to the sampling period of the system (i.e., τ 1 = t1 = T ); • the initial time-instants should be closer, while the successive ones can be sparser. Thus, given a certain number of time-instants c = p and the length of the horizon τ p in terms seconds, we set up a parabolic distribution of τ i and ti such that:

τ i = ti = τ p -T (p -1) 2 i 2 -2 τ p -T (p -1) 2 i + τ p -T (p -1) 2 + T (16)
for i = 1, . . . , p. In this way, the same predictive and control horizons can be reached with a reduced number of variables involved in the predictive model.

C. Satisfaction of the tasks

The satisfaction of the primary task along the predictive horizon can be expressed by imposing [START_REF] Chan | A weighted least-norm solution based scheme for avoiding joint limits for redundant joint manipulators[END_REF] at each time-instant of the horizon. This can be approximated by only considering the predictive timeinstants τ i and expressed as the minimization of the following cost function:

1 2 J * (q (τ) ) q(τ) -X (τ) s (τ) 2 (17)
where:

J * (q (τ) ) = blkdiag J q(t 0 + τ 1 ) , . . . , J q(t 0 + τ p ) ; q(τ) = q(t 0 + τ 1 ) T , . . . , q(t 0 + τ p ) T T ;

s (τ) = s(t 0 + τ 1 ), . . . , s(t 0 + τ p ) T ; X (τ) = blkdiag x (σ (t 0 ) + τ 1 ), . . . , x (σ (t 0 ) + τ p ) ;
The cost function [START_REF] Antonelli | A new on-line algorithm for inverse kinematics of robot manipulators ensuring path tracking capability under joint limits[END_REF] is nonlinear in the control action u (t) , since the future joint configuration q (τ) needed in the calculation of the Jacobians depends on u (t) . However, at each step, the prediction of the joint configuration can be calculated by means of ( 14) (with d = 0). Then, in the successive step, this prediction can be used to compute an estimation of the future Jacobians [START_REF] Tamas | The NEPSAC nonlinear predictive controller in a real life experiment[END_REF]. Thus, by considering such approximated Jacobians, ( 17) is now quadratic in q(τ) . By substituting [START_REF] Dahl | Torque-limited path following by online trajectory time scaling[END_REF] with d = 1, ( 17) is quadratic also in the control action u (t) . A quadratic cost function w 1 can therefore be defined as:

w 1 := 1 2 J * ( q)L 1 z 0 + J * ( q)F 1 u (t) -X (τ) s (τ) 2 ( 18
)
where q is the prediction of the joint configuration calculated at the previous time step by means of ( 14) (with d = 0). To recover from position errors, a CLIK strategy can be implemented as primary controller of a cascade architecture where the MPC is the secondary controller. In particular, the desired reference velocity X (τ) can be modified by adding a corrective term, that is, X (τ) = blkdiag x (σ (t 0 ) + τ 1 ) + K e c , x (σ (t 0 ) + τ 2 ) + K e c , . . . . . . , x (σ (t 0 ) + τ p ) + K e c ; [START_REF]Motion control of redundant robots under joints constraints: Saturation in the null space[END_REF] where K ∈ R m×m is a positive-definite diagonal matrix, and e c ∈ R m represents the current Cartesian position error, between the desired and the current positions. When the task orientation is expressed in terms of Euler angles and the analytical Jacobian is used in [START_REF] Flacco | Control of redundant robots under hard joint constraints: Saturation in the null space[END_REF], the Cartesian error can be simply defined as e c := x dx, where x d is the desired Cartesian position at the current time-instant, and x is the current Cartesian position. However, when angular velocities are used in the representation of the desired task velocity x , the Cartesian orientation error should be defined as the mutual orientation, in terms of axis-angle representation, between the desired and the current end-effector frames [START_REF] Siciliano | Six-degree-of-freedom impedance robot control[END_REF]. Further details about different orientation representations in the CLIK strategy can be found in [START_REF] Caccavale | The role of Euler parameters in robot control[END_REF].

Note that the CLIK controller bandwidth is set to be much smaller than the MPC controller one. This implies that the secondary loop can be reasonably approximated as a unitary gain. Thus, the stability of the primary loop can be guaranteed by imposing that K belongs to a region that depends on the sampling period [START_REF] Roveda | Highaccuracy robotized industrial assembly task control schema with force overshoots avoidance[END_REF].

The primary task [START_REF] Flacco | Control of redundant robots under hard joint constraints: Saturation in the null space[END_REF] can now be included in a priority framework with other secondary tasks, by following a hierarchical approach [START_REF] Kanoun | Kinematic control of redundant manipulators: generalizing the task priority framework to inequality tasks[END_REF]. The resulting hierarchical QP should (in order of priority): i) optimize the primary task (18); ii) limit the scaling; iii) optimize the lower-priority tasks. Considering a generic secondary task in the form E u (t) = g, this can be formalized as a lexicographic multi-objective quadratic problem, similarly to (9): lex min

u (t) , s (τ) w 1 ,w 2 , w 3 } subject to u (t) ∈ S , s i ∈ [0, 1] (20)
where w 2 := 1 2 1s (τ) 2 , w 3 := 1 2 E u (t)g 2 , and w 1 is defined in [START_REF] Flacco | Control of redundant robots under hard joint constraints: Saturation in the null space[END_REF]. The number of optimization variables involved in ( 20) is given by n c + p. In order to efficiently handle prioritized problems in the form of ( 20), fast dedicated QP solvers for hierarchical task have been recently developed, for instance, in [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF], while in [START_REF] Liu | Generalized hierarchical control[END_REF] a technique to solve hierarchical tasks in a single enlarged QP was proposed. Alternatively, a weighted approach can be adopted, in order to solve a single multi-objective QP with n c + p variables, in which the tasks have different weights, according to their priority. This can be written as:

u (t) s (τ) = argmin u∈S , s i ∈[0,1] λ 1 J * ( q) L 1 z 0 + J * ( q) F 1 u (t) -X (τ) s (τ) 2 +λ 2 1 -s (τ) 2 + λ 3 E u (t) -g 2 (21)
with λ 1 λ 2 λ 3 .

Examples of secondary tasks can be, for instance,

• joint range availability:

E = F 0 ; g = (Q max + Q min )/2 -L 0 z 0 ;
• minimum norm of the control action: E = I nc ; g = 0.

D. Analysis and design of the constraints

This section aims at translating the kinematic limits of the joints into linear constraints, in order to include them in the QP formulation. Note that, by means of the predictive equations ( 14), the bounds will be imposed along the whole control horizon. In this way, the optimal control action calculated at each step will take into account the presence of the constraints for the future time instants as well.

As introduced in Section II-B, the IK resolution should consider the boundaries of the kinematic variables, such as configuration, velocity and acceleration. Thus, by choosing the order r of the resolution algorithm, and by using [START_REF] Dahl | Torque-limited path following by online trajectory time scaling[END_REF], the following set of linear inequalities results:

Q (d) min -L d z 0 ≤ F d u (t) ≤ Q (d) max -L d z 0 ∀ d = 0, . . . , r -1 Q (r) min ≤ u (t) ≤ Q (r) max ( 22
)
The linear inequalities [START_REF] Hedjar | Nonlinear Receding-Horizon control of rigid link robot manipulators[END_REF] permit to take into account the joint boundaries up to the rth derivative of the configuration. However, the IK should consider the viability set condition [START_REF] Hollerbach | Dynamic scaling of manipulator trajectories[END_REF] as well, as mentioned in Section II-B. Unfortunately, ( 13) is nonlinear in the joint configuration and velocity and cannot be linearized as in the local resolution. In order to approximate it with a linear constraint a further analysis is required. First of all, it is important to point out that shrinking the state space (q, q) to a more conservative admissible set does not guarantee the satisfaction of the constraint in the future. Fig. 1: Admissible velocity area. Black: area given by the constraints ( 10)-( 13); Red: infeasible linearization of the viability condition; Blue: feasible linearization ( 23), [START_REF] Ferrara | A robust MPC/ISM hierarchical multi-loop control scheme for robot manipulators[END_REF] of the viability condition. Note that Proposition 1 yields to the graphical result ∆ 1 = ∆ 2 and ∆ 3 = ∆ 4 (if the boundaries are symmetric all the quantities are equal).

For instance, for every joint, consider the linear approximation of the constraint [START_REF] Hollerbach | Dynamic scaling of manipulator trajectories[END_REF] given by the line passing through the points (Q max,i , 0) and

(Q max,i + q2 i 2 Qmin,i
, Qmax,i ) (red line in Figure 1). When the constraint is activated, the state is forced to follow the trajectory given by the curve qi = 2 Qmin,i Qmax,i

q i -2 Qmin,i Qmax,i Q max,i . The corresponding ac- celeration is given by qi = d qi /dt = 2 Qmin,i Qmax,i qi , with qi ∈ [0, Qmax,i ]. Thus, max | qi | = max qi ∈[0, Qmax,i ] 2 Qmin,i Qmax,i qi = 2| Qmin,i |.
This means that the acceleration needed to move from the current value of qi to zero without violating the linear constraint is twice the maximum acceleration from which the constraint was constructed. Following a similar reasoning, it is possible to obtain the following proposition. Proposition 1: Consider the upper viability condition [START_REF] Hollerbach | Dynamic scaling of manipulator trajectories[END_REF] and its linear approximation in the form qi ≤ α q i + β . Imposing the passage of the curve through the point (Q max i , 0) of the (q i , qi ) plane, then the acceleration does not exceed the value Qmax,i if

α = - |ρ min,i | Qmax,i ; β = |ρ min,i | Qmax,i Q max,i (23)
where |ρ min,i | ≤ | Qmin,i |. Likewise, for the lower viability condition and its linear approximation qi ≥ γ q i + δ it results

γ = |ρ max,i | Qmin,i ; δ = - |ρ max,i | Qmin,i Q min,i (24)
where

|ρ max,i | ≤ | Qmax,i |.
Proof: Let the curve qi ≤ α q i + β be the approximation of the viability condition [START_REF] Hollerbach | Dynamic scaling of manipulator trajectories[END_REF]. Imposing the passage through the point

(Q max,i , 0) yields α Q max,i + β = 0 ⇒ α = -β /Q max,i . Therefore, qi = (-β / Qmax,i ) q i + β . The acceleration is given by qi = d qi /dt = (-β / Qmax,i ) qi , with qi ∈ [0, Qmax,i ]. Thus, max | qi | = max qi ∈[0, Qmax,i ] -β Q max,i qi = Qmax,i Q max,i | -β |.
Assuming β is negative and imposing max

| qi | ≤ | Qmin,i | gives β ≤ ( Qmin,i Q max,i )/ Qmin,i ⇒ α ≥ ( Qmin,i )/ Qmin,i or equivalently α = -|ρ min,i |/ Qmax,i and β = (|ρ min,i |/ Qmax,i ) Q max,i , where |ρ min,i | ≤ | Qmin,i |.
Once the constraint related to the viability set has been wisely linearized, its inclusion in the predictive framework is straightforward. As a matter of fact, it has to be:

qi (t) ≤ α i q i (t) + β i qi (t) ≥ γ i q i (t) + δ i ∀ t ∈ tc ∀ i = 1, . . . , n. (25)
Reminding the predictive equation ( 14), the constraints become:

(F 2 -ᾱF 1 ) u (t) ≤ ( ᾱL 1 -L 2 ) z 0 + β (F 2 -γ F 1 ) u (t) ≥ ( γ L 1 -L 2 ) z 0 + δ (26)
where

ᾱ := diag [α 1 , . . . , α n ], . . . , [α 1 , . . . , α n ] ∈ R n c×n c , β := [β 1 , . . . , β n ], . . . , [β 1 , . . . , β n ] T ∈ R n c
and likewise for γ and δ .

IV. NUMERICAL RESULTS

This section aims at testing the proposed method via simulation on a 9-dof redundant manipulator. The robot is composed of a 7dof KUKA iiwa14 manipulator, mounted on an actuated base, which consists in a 2-dof planar arm hanging on the ceiling. The manipulator is equipped with a 0.2m-long tool, aligned with the rotation axis of the last joint. The configuration and velocity limits are given by: The acceleration bounds are chosen equal to ± 12 rad/s 2 for all joints. The method is implemented with a sampling period T = 1 ms.

Q max = -Q min =[ 3.00,
Regarding the task, consider the two following reference paths:

• Path X1: a circle in the plane x = 0.35 m with center P 1 (0.35, -0.3, 1.3) m and radius R = 0.25 m. The circle has to be performed twice to complete the task; • Path X2: a sinusoidal curve z = 1.2 + 0.2 sin(4π -8πx) in the plane y = 0.1 m, from the initial point P 2 (0.5, 1, 1.2) m to the final point P 3 (0, 0.

m. For each of these reference paths, ten random paths are generated by perturbing the path parameters by adding a random number between ± 0.05 m (with uniform distribution).

The position of the end-effector should perform each path by following the nominal longitudinal timing law

η(σ ) = 6 T 5 tot σ 5 - 15 T 4 tot σ 4 + 10 T 3 tot σ 3 (27)
where T tot is the desired total duration of the nominal task. The following values of T tot are considered: T tot = 0.5, 1, 2, 4 s for Path X1, and T tot = 0.8, 1.2, 1.6, 2 s for Path X2. The longitudinal velocity profile of the end-effector can be obtained by differentiating [START_REF] Schuetz | Predictive online inverse kinematics for redundant manipulators[END_REF]. As secondary task, the robot elbow has to be kept as high as possible. This is expressed by defining a required velocity to the zcoordinate of the origin of the elbow frame in the form ẋ2 = |dη/dσ |. From a practical perspective, such behavior could be useful to reduce as much as possible the occupation of the space below the robot.

The following algorithms are compared:

• the local method (9) at acceleration level [START_REF]Motion control of redundant robots under joints constraints: Saturation in the null space[END_REF];

• the proposed predictive method, with number of control and predictive time instants c = p = 3 and predictive horizon τ p = 100 ms (namely, {τ i } = {1, 26, 100} ms); • the proposed predictive method, with c = p = 5 and predictive horizon τ p = 100 ms (namely, {τ i } = {1, 7, 26, 57, 100} ms). Referring to Section III-C, simulations related to Path X1 were performed by implementing the hierarchical approach, whereas simulations related to Path X2 were performed by implementing the weighted approach, with λ 1 = 10 6 , λ 2 = 10 3 , λ 3 = 10. The proportional gain related to the CLIK strategy is K = 200 I, where I is the identity matrix. Simulations were performed in Matlab.

The evaluated parameters are the deformation with respect to the desired path in terms of Euclidean distance with respect to the closest TABLE I: Numerical comparison between the local method and the proposed predictive method (p = 3 and p = 5). point on the desired path, the scaling factor along the task, and the satisfaction of the secondary task. The results are shown in Table I. For each method and for each value of T tot , the table shows the maximum error e max , the mean error e mean , the mean scaling factor s mean , and the average height of the robot elbow z elb,mean , given by the mean of the ten randomly-generated tasks. First of all, note that e max and e mean decrease for greater values of T tot (i.e., less demanding tasks). As regards the scaling values, as expected, more demanding tasks lead to smaller values of s mean , since a greater scaling of the timing law is needed. Consider now the comparison of the various methods. Firstly, e max and e mean given by the predictive methods always result to be significantly better than ones given by the corresponding local method. Secondly, the values of s mean given by the predictive methods are usually better than the ones given by the corresponding local method. An exception is represented by the fastest tasks. In such cases the predictive methods tend to perform a heavier scaling of the timing law to preserve the original geometrical path, which is the priority of the method. Simulations using Gazebo [START_REF] Koenig | Design and use paradigms for Gazebo, an open-source multi-robot simulator[END_REF] have then been performed to assess the effectiveness of the method in a realistic scenario. Gazebo receives the motor torque signals by a simple cascade P-PI control architecture with inertia matrix decoupling. The inner-and outer-loop bandwidths are respectively equal to 11 Hz and 4 Hz. In particular, few Cartesian tasks performed by the robot are shown in Figure 2. Table II shows the values of e max and e mean (considering also the controlled system tracking error), and s mean for such tests. In Figure 2a, it is possible to note that, for Path X1 and T tot = 2 s, the local method (dasheddotted blue line) gives a significant deformation of the desired path, while the predictive method (dashed red line) is able to avoid such deformation. This behavior is even more evident when the desired task becomes more demanding: as shown in Figure 2b for T tot = 0.5 s, the local method (dashed-dotted blue line) leads to a distorted path, especially in the final deceleration phase (zoomed window), while the predictive strategy (dashed red line) still gives satisfactory results. The same applies for Path X2, as shown in Figure 2c and2d.

As an illustrative example, the trends of the kinematic variables for Path X1, T tot = 2 s, and p = 3 are shown in Figure 3. First of all, the kinematic variables of the joint are always within their bounds. Moreover, as shown in the fourth column of the figure, the torque applied to the joints are kept within their limits, thanks to the appropriate choice of the acceleration limits.

Finally, the proposed method has also been tested at jerk level (r = 3). In many applications, for example to cope with unmodeled elasticities, it is advisable to reduce also higher-order derivatives to increase smoothness. By using r = 3, such soft constraints are taken into account as a part of the cost function. Joint accelerations are shown in Figure 4, which shows smoother trends of such variables, TABLE II: Numerical comparison between the local method and the proposed predictive method in Figure 2.

[s] [s] Local 8.3 • 10 -4 1.2 • 10 -4 Predictive (p = 3) 1.1 • 10 -3 2.1 • 10 -4 Predictive (p = 5) 2.3 • 10 -3 3.8 • 10 -4
since the algorithm also considers jerk minimization.

An important issue to take into account when using an optimization-based algorithm with fast sampling periods is the required computational time. For this reason, the mean and the maximum cycle time required by the algorithm are evaluated. Referring to the above-mentioned simulations, Table III shows the mean and the maximum computational time taken to perform Path X1. The simulations were performed by using the parametric active-set QP solver QpOASES [START_REF] Ferreau | qpOASES: A parametric active-set algorithm for quadratic programming[END_REF], [START_REF] Ferreau | An online active set strategy to overcome the limitations of explicit mpc[END_REF], on a standard laptop mounting a 2.5 GHz Intel Core i5-2520M processor. As expected, the computational time required by the predictive method is generally higher than the one required by the local method, due to the greater number of variables involved in the online optimization. Moreover, the computational time obviously grows with the number of control moves c = p. However, it is important to remind that the extension of the IK problem to the predictive framework originally led to a constrained nonlinear problem. The proposed method converts such problem into a QP problem, and the required computational time is still comparable with the one taken by the purely local method, and, thus, suitable for online implementation. Indeed, considering the example with p = 3, it is interesting to note that the maximum computational time is in the order of 1 ms. Furthermore, its mean value is around 0.21 ms, that is, less than the double of the local method, although the number of optimization variables is three times bigger. Finally, it is worth stressing that, in order to handle the real time execution of the control algorithm for which a fixed sampling period has to be selected, it is in any case possible to set a soft maximum CPU time, within which the optimization shall stop. From a practical perspective, the use of such setting does not give rise to a significant performance decay, if such maximum time is chosen properly, according to the problem to solve.

V. CONCLUSIONS

A new predictive task-scaling technique has been proposed in this paper. It allows solving the IK of a redundant manipulator taking into account the kinematic limits of the robot. The task-scaling strategy preserves the geometrical task when the desired motion is infeasible for the robot capabilities, dramatically reducing the geometrical error in case of very demanding tasks. Moreover, compared to local techniques, it is possible to obtain much better results in terms of task deformation and trajectory scaling. The use of a single-step linearization of the predicted task around the previous prediction of the joint configuration permits to convert the online optimization into a QP. In this way, the computational complexity is dramatically decreased, making the method suitable for online implementation.

17 •

 17 10 -1

Figure

  Figure e max [m] e mean [m] s mean Local p = 3 Local p = 3 Local p = 3 2a 1.66 • 10 -2 8.69 • 10 -4 1.11 • 10 -3 2.78 • 10 -4 0.909 0.918 2b 1.01 • 10 -1 7.33 • 10 -3 1.30 • 10 -2 4.21 • 10 -3 0.294 0.102 2c 3.09 • 10 -2 1.32 • 10 -3 2.79 • 10 -3 3.68 • 10 -4 0.855 0.854 2d 1.53 • 10 -1 5.31 • 10 -3 4.39 • 10 -2 5.84 • 10 -4 0.471 0.594

Fig. 2 :

 2 Fig. 2: Dotted black: ideal reference path. Dashed-dotted blue: path obtained by the local method. Dashed red: path obtained by the predictive method with r = 2, p = 3.

Fig. 3 :

 3 Fig.3: Solid black: joint configuration, velocity, acceleration and torque obtained by the predictive method with r = 2, p = 3, Path X1, T tot = 2 s. Dashed-dotted gray: joint configuration, velocity, acceleration and torque obtained by the local method, Path X1, T tot = 2 s. Dashed blue: joint configuration, velocity, acceleration and torque limits. The ninth joint is omitted for it does not affect the motion.

10 Fig. 4 :

 104 Fig. 4: Solid black: joint acceleration obtained by the predictive method with r = 3, p = 3, Path X1, T tot = 2 s. Dashed-dotted gray: joint acceleration obtained by the local method, Path X1, T tot = 2 s. Dashed blue: joint acceleration bounds. The ninth joint is omitted for it does not affect the motion.

  3.00, 2.96, 2.09, 2.96, 2.09, 2.96, 2.09, 3.05 ] T [rad] Qmax = -Qmin =[ 1.48, 1.48, 1.48, 1.48, 1.74, 1.30, 2.26, 2.35, 2.35 ] T [rad/s]

  1.01 • 10 -2 4.48 • 10 -3 1.33 • 10 -2 2.66 • 10 -3 2.22 • 10 -3 10 -2 2.14 • 10 -3 1.47 • 10 -3 2.17 • 10 -3 6.93 • 10 -4 5.85 • 10 -4 10 -2 2.59 • 10 -4 1.10 • 10 -4 7.22 • 10 -4 2.24 • 10 -5 1.97 • 10 -5

				0.286 0.101 0.082	1.28	1.59	1.65
	1	3.40 • 0.617 0.386 0.343	1.27	1.56	1.57
	2	1.67 • 0.920 0.964 0.962	1.33	1.39	1.40
	4	2.56 • 10 -5 2.40 • 10 -5 2.36 • 10 -5	5.59 • 10 -6 6.53 • 10 -6 6.34 • 10 -6	0.998 1.000 1.000	1.32	1.51	1.50
			Path X2
	0.8	1.52 • 10 -1 4.53 • 10 -3 1.50 • 10 -3	2.51 • 10 -2 2.15 • 10 -4 1.92 • 10 -4	0.479 0.556 0.455	1.30	1.39	1.38
	1.2	1.10 • 10 -1 1.51 • 10 -3 6.12 • 10 -4	1.29 • 10 -2 5.64 • 10 -5 6.02 • 10 -5	0.578 0.773 0.735	1.30	1.35	1.40
	1.6	3.43 • 10 -2 8.88 • 10 -4 4.05 • 10 -4	1.41 • 10 -3 2.20 • 10 -5 2.06 • 10 -5	0.869 0.928 0.909	1.31	1.41	1.37
	2.0	6.38 • 10 -4 2.76 • 10 -4 4.43 • 10 -5	1.29 • 10 -5 1.15 • 10 -5 1.19 • 10 -5	0.992 0.982 0.985	1.30	1.40	1.46

  8.69 • 10 -4 1.11 • 10 -3 2.78 • 10 -4 0.909 0.918 2b 1.01 • 10 -1 7.33 • 10 -3 1.30 • 10 -2 4.21 • 10 -3 0.294 0.102 2c 3.09 • 10 -2 1.32 • 10 -3 2.79 • 10 -3 3.68 • 10 -4 0.855 0.854 2d 1.53 • 10 -1 5.31 • 10 -3 4.39 • 10 -2 5.84 • 10 -4 0.471 0.594

TABLE III :

 III Mean and maximum computational time taken by the local and the predictive (p = 3 and p = 5) methods for Path X1.

	Method	t max	t mean

ACKNOWLEDGMENT

The research leading to these results has been partially funded by the European Union H2020-CS2-738039-Eureca: Enhanced Human Robot cooperation in Cabin Assembly tasks.


Mario Mekhael

Aren Chris Labaki

Joe Bizdikian

Ziad Bakouny

Fares Yared Joeffroy Otayek

Fares Yared a , Abir Massaad

Wafa Skalli
email: wafa.skalli@ensam.eu

Ismat Ghanem

Ayman Assi
email: ayman.assi@usj.edu.lb

Aren Joe Bizdikian

How do skeletal and postural parameters contribute to maintain balance during walking?

Keywords:

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Balance maintenance depends on multiple systems such as the vestibular, proprioceptive, visual and musculoskeletal systems [START_REF] Maurer | Multisensory control of human upright stance[END_REF]. The concept of the conus of economy presents the variability in standing posture which allows the maintenance of static equilibrium while minimizing energy expenditure and fatigue on the muscles of the back and the lower limbs. An unstable standing posture can lead to subsequent pain and disability (Jean Dubousset, 1994). In dynamic, balance consists in minimizing energy expenditure and avoiding fall during walking [START_REF] Winter | Human blance and posture control during standing and walking[END_REF].

Several techniques for evaluating balance during gait have already been proposed. This has been evaluated either by determining local and orbital dynamic stability [START_REF] Dingwell | Differences between local and orbital dynamic stability during human walking[END_REF]) by using tri-axial accelerometry [START_REF] Auvinet | Reference data for normal subjects obtained with an accelerometric device[END_REF] or by studying the center of mass/center of pressure COM-COP relationship [START_REF] Lee | Detection of gait instability using the center of mass and center of pressure inclination angles[END_REF]. The center of mass (COM) is defined as the barycenter of the subject, while the center of pressure (COP) is the application point of the ground reaction force captured by forceplates during gait. The angle between the COM-COP and the vertical can be evaluated in both frontal and sagittal planes, and its increase is correlated with increased instability. During the gait cycle, the sagittal COM-COP angle presents large variability in contrast with the frontal COM-COP, whose variations are limited [START_REF] Lee | Detection of gait instability using the center of mass and center of pressure inclination angles[END_REF]. The variation of these two angles during gait has been shown to reflect dynamic gait balance [START_REF] Lee | Detection of gait instability using the center of mass and center of pressure inclination angles[END_REF][START_REF] Paul | Gait stability improvement after fusion surgery for adolescent idiopathic scoliosis is influenced by corrective measures in coronal and sagittal planes[END_REF].

While demographic parameters (age, sex, weight and height) are known to affect gait [START_REF] Gao | Dynamic walking stability of elderly people with various BMIs[END_REF][START_REF] Hamacher | Gait stability and its influencing factors in older adults[END_REF], their relationship with gait stability is still unelucidated. Moreover, previous studies have determined that skeletal and postural alterations can affect gait balance. For instance, the correction of coronal and sagittal spinal malalignment in scoliotic patients decreased mediolateral excursion during gait [START_REF] Paul | Gait stability improvement after fusion surgery for adolescent idiopathic scoliosis is influenced by corrective measures in coronal and sagittal planes[END_REF]. In another setting, patients with adult spinal deformity (ASD) are known to have alteration of their postural malalignment that affects not only their spine and pelvis but also the positioning of their head and lower limbs; a chain of compensatory mechanisms is established in order to maintain a balanced standing posture in the conus of economy [START_REF] Diebo | Adult spinal deformity[END_REF]. Patients with ASD are known to have altered quality of life that affects their daily life activities such as walking. Several studies started to highlight gait abnormalities in patients with ASD (Assi et al., 2019;[START_REF] Haddas | The use of gait analysis in the assessment of patients a ffl icted with spinal disorders[END_REF]. However, it is still unknown how their skeletal and postural deformities are related to their gait balance. Any investigation on the relationship between skeletal or postural deformities and gait balance in the setting of ASD or other musculoskeletal pathologies would benefit from a baseline or normative relationship that is still unknown.

As a first step, the aim of this study was to evaluate the normative relationship between gait balance and both demographic and skeletal and postural parameters in asymptomatic subjects.

Materials and methods

Study design

This is a cross-sectional IRB approved study (CEHDF285-2016) evaluating the different determinants of balance during gait in asymptomatic adult subjects recruited from our institution. All subjects signed an informed consent form. Subjects were excluded if they presented with prior history of self-reported back pain or orthopedic surgery. Exclusion criteria included also the presence of scoliosis (frontal Cobb > 10°), Scheuermann's kyphosis, leg length discrepancy or any radiologic determinant of adult spinal deformity [START_REF] Pellisé | Impact on health related quality of life of adult spinal deformity (ASD) compared with other chronic conditions[END_REF] detected on radiographs.

Data collection

Demographics (sex, age, weight and height) were collected. Subjects underwent three-dimensional gait analysis (3DGA) using a Vicon (Vicon Motion Systems, Oxford, UK) optoelectronic motion system (7 MX3 infrared cameras, 200 Hz). Marker placement was based on the modified Helen Hayes protocol and was applied as recommended in the Plug-in-Gait® model. Subjects were asked to walk along a 10-m walkway, at a self-selected speed, equipped with two calibrated forceplates (OR6-7-1000 AMTI®, MA, USA). Several trials were recorded and only those where the subjects placed a full foot on each of the forceplates were considered. The pipeline in Workstation® (Vicon Motion Systems, Oxford, UK) was used to process data: fill gap routine ± 10 frames and Woltring filter with a scale of 10. Consistent kinematic and kinetic curves were validated using Polygon® (Vicon Motion Systems, Oxford, UK) and inconsistent trials were eliminated. One representative trial was then considered.

The center of mass (COM) was determined as the vector sum of the individual segment mass vectors, estimated from the kinematic centroid obtained from the Plug-In-Gait model. The center of pressure (COP) was determined as the application point of the ground reaction force obtained from the forceplates. The COP data were time-synchronized with motion data. During the double-stance phase, when the COP changed between the forceplates, a resultant COP was calculated for both feet [START_REF] Jian | Trajectory of the body COG and COP during initiation and termination of gait[END_REF].The COM-COP angle with the vertical was calculated in both the frontal and sagittal planes during the gait cycle for each subject [START_REF] Lee | Detection of gait instability using the center of mass and center of pressure inclination angles[END_REF]) (Fig. 1).

In the frontal plane, the maximum of the COM-COP angle was chosen as the highest value for the left or the right side, depending on which was higher. As for the sagittal plane, the maximum value was considered to be the highest positive value and therefore when the COP was in front of the COM. The minimum of the COM-COP was the lowest negative value, therefore the COP being behind the COM.

The maxima, minima and range of motion (ROM) of the COM-COP angle in the frontal and sagittal planes were considered as parameters of gait balance and were calculated over the whole gait cycle for each subject. An increase in any of these parameters would indicate a decrease in balance.

Following the gait analysis, subjects underwent a full body biplanar X-ray exam (EOS Imaging®, Paris, France). All subjects were placed in the consensual free-standing position with shoulders flexed at 45°and hands placed on the zygomatic processes, and with slightly shifted feet [START_REF] Chaibi | Fast 3D reconstruction of the lower limb using a parametric model and statistical inferences and clinical measurements calculation from biplanar X-rays[END_REF]. This position was adopted to avoid the superimposition of subjects' arms over their spines, or the left and right femoral condyles, on lateral radiographs. Three-dimensional reconstructions of the spine and the lower limbs were obtained from the biplanar X-rays using SterEOS® (EOS imaging, Paris, France). Three-dimensional reconstruction of the pelvis was obtained using a specific software (Arts et Métiers ParisTech, Paris, France). Then, skeletal and postural radiological parameters of the spine, pelvis, lower limbs and global posture were calculated: sagittal vertical axis SVA (center of C7 vertebral body plumbline to the posterior corner of the sacral plate), CAM plumb line (center of auditory meatus plumbline to the middle of the hip axis), thoracic kyphosis (T1T12), lumbar lordosis (L1S1), pelvic incidence, pelvic tilt, sacral slope [START_REF] Bendaya | Healthy vs. osteoarthritic hips: A comparison of hip, pelvis and femoral parameters and relationships using the EOS® system[END_REF][START_REF] Dubousset | Skeletal and spinal imaging with EOS system[END_REF][START_REF] Dubousset | EOS stereo-radiography system: Whole-body simultaneous anteroposterior and lateral radiographs with very low radiation dose[END_REF][START_REF] Duval-Beaupere | Sagittal profile of the spine prominent part of the pelvis[END_REF][START_REF] Legaye | Pelvic incidence: A fundamental pelvic parameter for three-dimensional regulation of spinal sagittal curves[END_REF]Jean Legaye & Duval-Beaupère, 2005;[START_REF] Vialle | Radiographic analysis of the sagittal alignment and balance of the spine[END_REF], femoral neck length, neck shaft angle, femoral torsion, valgus/varus of the knee, femoral anteversion, tibial torsion, tibial mechanical angle, acetabular abduction, acetabular anteversion, acetabular tilt, vertical center edge angle (VCE) [START_REF] Anda | The acetabular sector angle of the adult hip determined by computed tomography[END_REF][START_REF] Stem | Computed tomography analysis of acetabular anteversion and abduction[END_REF][START_REF] Zilber | Variations of caudal, central, and cranial acetabular anteversion according to the tilt of the pelvis[END_REF], anterior (AASA) and posterior (PASA) acetabular sector angles [START_REF] Anda | The acetabular sector angle of the adult hip determined by computed tomography[END_REF], percentage of acetabular coverage by the femoral head [START_REF] Humbert | 3D evaluation of the acetabular coverage assessed by biplanar X-rays or single anteroposterior Xray compared with CT-scan[END_REF] and the sacro-acetabular angle (Fig. 2).

Statistical analysis

Descriptive statistics (mean and standard deviation) were calculated for the COM-COP angles in both the sagittal and frontal planes, as well as for the spino-pelvic, hip and lower limb parameters.

In order to evaluate the relationship between the gait balance (COM-COP) variables, and demographic and skeletal parameters, a univariate analysis was conducted using Pearson's correlation test. The level of significance was set at 0.05. P-values were adjusted using Benjamini-Hochberg procedure to account for multiple testing (n = 312 correlations).

Then, in order to evaluate the main determinants of gait balance from demographic and skeletal parameters while accounting for correlations between predictor variables, a multivariate analysis was conducted using a stepwise multiple linear regression (SMLR). Probability for entry was 0.05 / Probability for removal was 0.1. The dependent variables were the COM-COP angle parameters in the sagittal and frontal planes; the independent variables were demographics and skeletal parameters that showed significant correlations to COM-COP variables in the univariate analysis, as well as other non-significant parameters that were judged to be clinically important to be included in the model. In total, 5 models were run for the following COM-COP parameters: maximum frontal angle, ROM frontal angle, minimum sagittal angle, maximum sagittal angle and ROM sagittal angle. The level of significance of the model was adjusted using Bonferroni correction to account for multiple models (n = 5; thus, level of significance was 0.01 instead of 0.05). Moreover, the level of significance for each predictor variable was adjusted using the Benjamini-Hochberg procedure in order to account for the number of predictor variables included in the model (n = 18). The validity of multiple linear regression assumptions was checked. The adjusted R 2 of the models as well as β and p-values of predictor variables were reported.

Statistics were undertaken using Xlstat® (version 2019.1.2, Addinsoft, Paris, France).

Results

In total, 115 asymptomatic subjects (age: 29 ± 7 years [18-59], weight: 70 ± 9 Kgs, height: 169 ± 10 cm, 57F: 58 M) were included in this study.

Descriptive statistics

Subjects had an average pelvic tilt of 12°with a pelvic incidence of 48°. Their L1S1 lumbar lordosis had a mean of 60°associated with a T1T12 thoracic kyphosis of 45°. Their C7 was positioned in the sagittal plane behind the sacrum (SVA = -12 cm) as well as their center of auditory meatus (CAM-HA = -25 cm). They had a pelvic width of 25 cm on average. The acetabulum was tilted with an average of 23°. The abduction had a mean of 55°and the anteversion was 17°. The acetabulum was covered posteriorly with a mean of 96°(PASA) and anteriorly 60°(AASA). The lateral coverage of the acetabulum (VCE) had a mean value of 30°(detailed values and other parameters are displayed in Fig. 3).

In the frontal plane, the maximum of the COM-COP reached 4.6 ± 1.3°and the ROM was 9.0 ± 2.4°. In the sagittal plane, the maximum of the COM-COP angle reached 19.7 ± 2.8°, the ROM was 34.1 ± 4.8°and the minimum was -14.4 ± 3.21°( Fig. 4).

Univariate analysis

As a first step in assessing the relationship between demographics and skeletal parameters with gait balance, simple correlations were computed between the variables.

The COM-COP angles were significantly correlated to the following demographic and skeletal parameters (Fig. 5): An increase of weight (r = 0.53, Fig. 6), age (r = 0.22 to 0.28), and height (r = 0.35 to 0.37) was correlated to an increase in the maximum and ROM of the frontal COM-COP angle. Regarding skeletal and postural parameters, an increase of a forward back (SVA: r = 0.2 to 0.23), thoracic kyphosis (T1T12: r = 0.22 to 0.24) and pelvic width (r = 0.24 to 0.25) was correlated to an increase in the maximum and ROM of the frontal COM-COP.

However, in the sagittal plane, an increase of height (r = -0.28 to -0.33, Fig. 6) was correlated to a decrease of maximum and ROM of the sagittal COM-COP. An increase of acetabular tilt (r = 0.25) and acetabular anteversion (r = 0.21) were correlated to an increase of the maximum sagittal COM-COP. In addition, a decrease in anterior acetabular coverage (r = -0.21) was correlated to an increase of the ROM sagittal COM-COP.

Multivariate analysis

Furthermore, in order to better investigate the main determinants of gait balance parameters among demographic and skeletal and postural multiple linear regressions were computed. The following parameters that were shown to be significant in the univariate analysis as well other parameters that were judged to be clinically significant were included in the multivariate models as dependent variables: age, weight, height, sex, SVA, CAM-HA, L1S1, T1T12, pelvic tilt, pelvic incidence, VCE, percentage of acetabular coverage over the femoral head, sacro-acetabular angle, acetabular tilt, acetabular abduction, acetabular anteversion, PASA, AASA.

Instability during gait was determined by the increase of the following parameters in each plane as shown in the multivariate analysis:

In the frontal plane, an increase in the maximum of the COM-COP angle was determined (p < .001) at 39% (adjusted R 2 ) by an increasing weight (β = 0.48, p < .001), being a male (β = 0.42, p < .001) and an increasing posterior acetabular coverage (β = 0.22, p = .006). An increasing ROM of the frontal COM-COP angle was determined (p < .001) at 35% by an increasing weight (β = 0.32, p = .002), being a male (β = 0.37, p = .001), and an increasing posterior acetabular coverage (β = 0.22, p = .007). In the sagittal plane, an increasing maximum of the COM-COP angle was determined (p < .001) 13% by a decreasing height (β = -0.38, p < .001) and an increasing SVA (β = 0.19, p = .004). An increasing ROM of the sagittal COM-COP was solely determined (p = .01) at 7.3% by a decreasing height (β = -0.28, p = .002).

Discussion

In this study, the examination of skeletal parameters and gait of 115 asymptomatic adult subjects showed that gait balance in the frontal and sagittal planes was influenced by both anthropometric and skeletal parameters. Frontal imbalance appeared to be more pronounced with increasing weight, height, age, SVA, T1T12 and pelvic width. Sagittal imbalance was found to be associated with an increased acetabular tilt, acetabular anteversion and a decreased anterior acetabular coverage. In addition, sagittal COM-COP was shown to increase with decreasing height and weight.

The hip skeletal parameters of our population were found to be similar to those previously reported in studies concerning asymptomatic adults [START_REF] Tannast | What are the radiographic reference values for acetabular under-and Overcoverage[END_REF]. Similarly, spino-pelvic parameters of this study were closely comparable with the results shown in previous studies [START_REF] Bakouny | Normative spino-pelvic sagittal alignment of Lebanese asymptomatic adults: Comparisons with different ethnicities[END_REF].

The COM-COP parameters of our population were also similar to the ones reported in previous studies [START_REF] Lee | Detection of gait instability using the center of mass and center of pressure inclination angles[END_REF].

The female sex appeared as being an important determining factor (adjusted R 2 > 0.33) of both the maximum (β: -0.426, M being the reference) and the ROM (β: -0.371) of the frontal COM-COP. It is known that the female pelvic anatomy is significantly different compared to males [START_REF] Kersnic | Determination of the femoral and pelvic geometrical parameters that are important for the hip joint contact stress: Differences between female and male[END_REF][START_REF] Wang | Gender differences in hip anatomy: Possible implications for injury tolerance in frontal collisions[END_REF]: the female pelvis being larger and wider than the male pelvis which is taller, narrower, and more compact. This anatomic property could influence balance during gait, thereby possibly explaining the increased frontal instability reported in women.

An increased weight appeared to be related to frontal imbalance. It was correlated to both the maximum (r = 0.53) and the ROM of the frontal COM-COP (r = 0.53). These results were reinforced by the multivariate analysis which showed that weight is a determinant of both frontal maximum and ROM COM-COP angles. These findings are further illustrated by a comparison of two subjects with differing weights who present important variations in their maximum frontal COM-COP angles (Fig. 7).

On the other hand, pelvic width was also correlated to both the maximum (r = 0.25) and the ROM (r = 0.24) of the frontal COM-COP angle suggesting that a wider pelvis morphology may be associated with an increased imbalance in the frontal plane.

Furthermore, the sagittal vertical axis is a commonly used parameter to study global alignment. The subject is considered to have sagittal malalignment if the SVA drops in front of the femoral heads [START_REF] Jackson | Radiographic analysis of sagittal plane alignment and balance in standing volunteers and patients with low Back pain matched for age, sex, and size: A prospective controlled clinical study[END_REF]. This parameter was used in the several classifications for adult spinal deformity subjects as an inclusion criteria for subjects having an SVA of over 4 or 5 cm [START_REF] Pellisé | Impact on health related quality of life of adult spinal deformity (ASD) compared with other chronic conditions[END_REF][START_REF] Schwab | Scoliosis Research Society-Schwab adult spinal deformity classification: A validation study[END_REF]. Glassman showed that an increased forward bending of the trunk is associated with increased pain and decreased function [START_REF] Glassman | The impact of positive sagittal balance in adult spinal deformity[END_REF]. The magnitude of these symptoms was also found to increase with augmented imbalance, i.e. when the SVA increases. An increased SVA is well known to affect a subject's quality of life and is therefore associated with altered HRQOL measures [START_REF] Kim | The cranial sagittal vertical axis (CrSVA) is a better radiographic measure to predict clinical outcomes in adult spinal deformity surgery than the C7 SVA: A monocentric study[END_REF]. In this study, the SVA was found to be positively correlated with the maximum (0.23) and the ROM (0.20) of the frontal COM-COP angles. SVA was also shown to be a determinant of maximum sagittal COM-COP angle (β = 0.19). Thus, an increase of SVA, that is usually seen in elderly people and especially in adults with spinal deformity [START_REF] Diebo | Adult spinal deformity[END_REF], is related to increasing imbalance during gait in both frontal and sagittal planes (Assi, Kawkabani, et al., 2019). Future studies should investigate about causes of gait imbalance in subjects with adult spinal deformity (ASD) and to verify if surgical correction of SVA can improve gait stability in these patients.

In addition to global posture, several acetabular parameters were found to affect gait, and especially in the sagittal plane. An increased acetabular tilt (r = 0.25) and acetabular anteversion (r = 0.21) were found to be correlated with increased maximum Fig. 6. Scatter plots of selected correlations between COM-COP angles and demographic and skeletal parameters.

sagittal COM-COP angle. Also, a decreased acetabular anterior coverage (r = -0.21) is correlated to an increased sagittal ROM of the COM-COP angle. These findings are further by a comparison of two subjects with considerable variations in their acetabular parameters who present with different maximum sagittal COM-COP angles (Fig. 8). Additionally, posterior acetabular coverage appeared to be an important determinant (adjusted R 2 > 0.35) of both maximum (β = 0.22) and ROM (β = 0.22) of the frontal COM-COP angle. These results show that when the acetabulum is more tilted in the sagittal plane, thus decreasing acetabular anterior coverage and increasing posterior acetabular coverage, this could contribute to increase sagittal imbalance during gait. These modifications in acetabular orientation are seen in subjects with adult spinal deformity due to their increased pelvic tilt that contributes to an increase of their acetabular tilt associated with a decreased anterior coverage and an increase of the posterior coverage of the (Assi et al., 2019). Thus, these acetabular orientation parameters might be related to the gait imbalance in subjects with adult spinal deformity. Future studies should test this hypothesis in patients with ASD.

Moreover, these results suggest that the orientation of the proximal hip, which forms the anatomical hinge between the head and trunk segments, where the centre of mass is located, and the lower limbs, which are the principal actors in gait, is a very important factor that correlates with sagittal dynamic imbalance, and therefore instability during gait. This is the first study to evaluate the normative relationship between demographic and skeletal parameters with gait stability. This study has some limitations; it did not take into account the muscular component when assessing the determinants of gait stability. Future studies should include the muscular factor, using EMG or musculoskeletal models, when evaluating stability during walking. As discussed earlier, future studies should also investigate gait balance in subjects with skeletal deformities in the spine and the hip such as in subjects with adult spinal deformity.

In conclusion, this study evaluated the demographic and skeletal determinants of gait balance in asymptomatic adult subjects. The results suggest that in addition to demographic parameters, acetabular orientation and SVA are important determinants of both frontal and sagittal balance during gait. Future studies should analyze the stability during gait in patients with spinal or hip deformities, therefore allowing an understanding of the deformity and its repercussion on walking, leading to better patient-oriented treatment options.

Fig. 1 .

 1 Fig. 1. Representation of COM-COP angle calculation in the sagittal and frontal planes.

Fig. 2 .

 2 Fig. 2. Skeletal and postural parameters calculated from 3D reconstructions obtained from full-body biplanar radiographs.

Fig. 3 .

 3 Fig. 3. Spino-pelvic, hip and global postural parameters calculated in 115 asymptomatic adults.

Fig. 4 .

 4 Fig. 4. Sagittal and frontal COM-COP angles calculated in 115 asymptomatic subjects.

Fig. 5 .

 5 Fig. 5. Correlation diagram between COM-COP angles and skeletal parameters in 115 asymptomatic subjects.

Fig. 7 .

 7 Fig. 7. Example of the maximum frontal COM-COP angle in two subjects with different weights.

Acknowledgment

This research was funded by the University of Saint-Joseph (grant FM183). The funding source did not intervene in study design; in the collection, analysis and interpretation of data; in the writing of the report; and in the decision to submit the article for publication.

Declarations of Competing Interest

None.

Author Statement

All the authors certify that they have seen and approved the final version of the manuscript being submitted. They warrant that the article is the authors' original work, hasn't received prior publication and isn't under consideration for publication elsewhere.