# Weight Lifting Trajectory Optimization for Variable Stiffness Actuated Robot

Authors: Chen Ji; Minxiu Kong; Ruifeng Li
Journal: Complexity
 (2019)

---

## Abstract

Variable stiffness actuator is a new actuation of robot which is inspired by motor control of human arm. And it is promising way to exploit the human-like performance and human-like motion. However, due to the mechanical complexity and redundancy of the actuators, it is not trivial to control the variable stiffness actuated robot to perform a human-like motion. In this paper, the weight lifting problem with variable stiffness actuated robot is studied. Then, the original problem is formulated as a constrained optimal control problem and transformed as a general nonlinear programming problem with Gauss pseudospectral method. Simulation studies demonstrate the effectiveness of the proposed method compared with the iLQR approach. Furthermore, the simulations results are presented to show the influence of stiffness variation of variable stiffness actuators on the weight lifting problem.


---

## Body

## 1. Introduction

The research on how to make the robots perform and move like humans has been paid great attention to in the past decades. It is a big challenge for robotic researchers to achieve human-like motion and performance, including but not limited to the design of biological mechanism and control method with the diverse requirement of performance for different scenarios. Inspired by motor control of human, a set of biomimetic actuations called variable stiffness actuators has been concentrated on [1, 2]. The main concept of these actuators is to modulate the output stiffness mechanically, instead of the software-based control method, like impedance control method [3]. Therefore, the redundant motors are introduced in the actuators in which the output position and stiffness can be modulated simultaneously. The introduced double actuation in the variable stiffness actuators can perform beyond the traditional actuation, like in the case of explosive movement [4, 5] and safe physical human-robot interaction. Besides, the potential applications of variable stiffness actuators to achieve human-like motion and performance still should be studied and extended. In this paper, weight lifting problem with variable stiffness actuated robot is focused on. Weight lifting in Olympic Games is the athlete attempts to lift weight, as heavier as possible, from the ground to the top of the head. In this paper, the weight lifting problem of robot is presented as payload maximization along a point-to-point motion in the minimum time and energy cost under dynamic and kinematic constraints in Section 3.In the literature, there are amounts of researches on payload maximization of different robots, such as industrial robot [6], humanoid robot [7], mobile manipulator [8], and compliant robot [9]. Nevertheless, the optimal control method is still the most common tool to solve this kind of problem, including indirect method based on Pontryagin’s maximum principle and direct parameter optimization [6]. Besides, some other methods are also employed to solve the weight lifting problem of robot, such as trial-and-error learning method [10], redundancy resolution method [11, 12], and polynomial-based trajectory optimization [13]. Thus, the optimal control method is a promising method to solve the weight lifting problem of variable stiffness actuated robot. However, the redundancy and strongly nonlinearity of variable stiffness actuators complicate the problem greatly. Considering the complexity of the system plants, it is very necessary to choose a suitable optimal control method.For the solutions of nonlinear constrained optimal control problem, there are plenty of researches on the optimal control method which can be categorized as indirect method, direct method, and dynamic programming method. And also some advanced optimization method, such as multiobjective optimal control method [14–16] and genetic algorithm, can also be applied for this problem. The main principle of indirect approach is to reconstruct the original optimal control problem as a two-point boundary value problem (TPBVP) based on Pontryagin’s maximum principle and first-order optimality conditions. Then, the TPBVP can be solved with traditional shooting method, such as single shooting method and multiple shooting method [17]. The main advantage of indirect method is its generality for the optimal control problem. However, it is suffered with difficult initialization of control and state variables and nonderivative plant dynamics. The successive approximation method which is one branch of dynamic programming is very popular for the complex optimal control problem for its computational efficiency. The iterative linear quadratic regulator/Gaussian (iLQR/G) approach [18] is utilized for the optimal control problem of variable stiffness actuated robot [4, 5]. Under the iLQR/iLQG framework, the original optimal control problem is discretized as the suboptimal problem which can be iteratively solved to obtain the optimal control law. However, solving an optimal control problem with state inequality constraints and equality constraints with the iLQR/iLQG approach is often computationally demanding and complicates the numerical treatment [19]. An alternative method is direct method which employs control and/or state parametrization to transform the original optimal control problem to a nonlinear programming problem (NLP) [20]. Instead of finding a solution to the first-order optimality conditions, the original infinite dimensional constrained functional minimization is directly converted to a finite dimensional function minimization that is easier to solve numerically. Based on different parameterization method, the direct method can be categorized into two main types, namely, direct collocation method and spectral method. In the direct collocation method, piecewise polynomials can be used to approximate the system dynamics at collocation points, like Euler method, Runge-Kutta method, and Trapezoid method. In the spectral method, the control and states can be parameterized using global polynomials, like Gauss pseudospectral method and Legendre pseudospectral method. However, Gauss pseudospectral method is preferred due to its high precision and fast rate of convergence [21]. Therefore, Gauss pseudospectral method becomes a very effective way in solving complex optimal control problems [22, 23] and is widely used in solving the problem of aircraft optimization [24].In this paper, Gauss pseudospectral method is used to solve the weight lifting problem for variable stiffness actuated robot and is presented in detail in Section4. Compared with the iLQR approach in the authors’ prior work [19], the simulation results on the same case show the proposed method can generate more precise and smoother trajectory. To the authors’ best knowledge, it is the first time that the weight lifting problem in case of variable stiffness actuated robot is reported. And the contribution of this work is the application of the Gauss pseudospectral method for the case of variable stiffness actuated robot. The proposed method is a practical method and can be used by the engineers or the researchers to implement trajectory planning efficiently.The remainder of this paper is as follows. The model of variable stiffness actuated robot with dynamic payload is presented in Section2. In addition, a roller-cam-based design of variable stiffness actuator is introduced in Section 5. To validate the method, simulations are performed for the case study of the two-link variable stiffness actuated robot.


## 2. Dynamic Model with Dynamic Payload of Variable Stiffness Actuated Robot

### 2.1. Link-Side Dynamic of Robot

The variable stiffness actuated robot can be treated as complaint joints with stiffness variation. The general dynamic model of rigid robot with compliant joints and dynamic payload is recalled here. Without loss of generality, consider the variable stiffness actuated robot as an open kinematic chain havingN+1 rigid links, interconnected by N joints, and actuated by M electrical motors. The link-side coordinate q∈RN describes the joint angles and the motor-side coordinate θ∈RM denotes the electric motor angles. Obviously, the variable stiffness actuated robot is redundantly driven by electrical motors. These motors can be named as position motor and stiffness motor which are in charge of the position modulation and for output stiffness modulation of robot link, respectively. Therefore, the motor-side coordinate is defined as θ=[θpos,θstiff]. And θpos∈RN are the motor angles for position motor and θstiff∈RM-N are the motor angles of stiffness motor.Then, the kinetic energy, gravitational potential energy, and elastic potential energy are defined using the same formalism in [25]. Following the Lagrangian approach, the equation of motion is derived and represented as follows.( 1 ) Kinetic energy: the total kinetic energy of the rigid body with dynamic payload is given by(1)Ek=∑i=1NKi+12mLv0L22,where dynamic payload is denoted asmL and is considered as lumped mass which is fixed at the end-effector of robot. v0L is the velocity vector of the payload which is relative to the base coordinate. Ki represents the kinetic energy of each rigid link.( 2 ) Gravitational potential energy: The gravitational potential energy depends on the configuration of the rigid body dynamics. And the center of the mass of the actuators is assumed that will not change due to the self-rotation of the rotors and the gears.(2)UG=∑i=0N-mig·r0ci-mLg·r0cL,wheremi is the mass of robot’s i-th rigid link and r0ci is the mass center vector of i-th rigid link relative to the base coordinate. r0cL is the mass center vector relative to the base coordinate of dynamic payload which is defined above.( 3 ) Elastic potential energy: Due to the elastic element in variable stiffness actuators, the conservative potential energy can be represented as follows:(3)UE=∑i=0N∫τiθstiff,φθ,qdφ,whereτ(θstiff,φ)∈RN=r(θstiff)F(φ,θstiff) represents the general form of the elastic torque vector. r∈RN×p denotes the position matrix of the moment arm and F∈Rp denotes the corresponding forces due to the elastic elements (normally is the spring). Based on Spong’s model of flexible joint of robot, the variable stiffness actuator could be considered as complaint actuator with virtual torsional spring and then the deformation of spring is represented as φ(θ,q) and simplified as φ in the rest of the paper. In general, there are two different ways to modulate the output stiffness of the variable stiffness actuator: one way is through the adjustment of the spring preload [2] and the other way is through modulating the transmission ratio between input and output torque, such as lever principle [26].Following the Lagrangian approach, the equation of motion is derived and represented as follows:(4)L=EK-UG-UE(5)ddt∂L∂q˙-∂L∂q=0.Then, the dynamic model of variable stiffness actuated robot with dynamic payload is represented as(6)Mq,mLq¨+Cq,q˙,mLq˙+Dqq˙+Gq,mL=τEφ,θstiff,whereM(q,mL)∈RN×N is the symmetric and positive definite inertia matrix of the rigid body dynamics. C(q,q˙,mL)∈RN×N is the matrix that represents the Coriolis term, Dqq˙∈RN are the forces due to viscous friction, G(q,mL)∈RN are the gravitational forces and influenced only by the configuration of the robot, and τE(φ,θstiff)∈RN are the elastic joint torques that affect the rigid body dynamics.


### 2.2. Motor-Side Dynamic Model of Robot

The redundant motors are responsible for the output position and stiffness of robot links. Through the control of motor angles, we can obtain the stiffness and position of robot links simultaneously. Considering the motor-side dynamic model which is coupled with link-side dynamics will complicate the solution of the optimal control problem. Fortunately, for the servomotors with a high transmission ratio or high-gain feedback position controllers [27], the motor-side dynamics can be considered as decoupled from the link side and are expressed as follows:(7)JMθ¨+Dθθ˙=τEφ,θstiff+τM,whereJM∈RM×M represents the inertia of the motors which is a constant and diagonal matrix. Dθθ˙∈RM represents the viscous friction forces in the drive train. τM∈RM is the output torque of servomotors. Then, the classical second-order position controller of servomotors can be designed as(8)τM=τ^q,θ-ξθ˙-κθ-θd,whereθd∈RM is the reference position of motors for the control loop. ξ=diag⁡{2αi}∈RM×M, κ=diag⁡{αi2}∈RM×M are user defined gains as critically damped controller which leads to a fast response without overshooting control effort. And τ^=-τ^E(φ,θstiff)+D^θθ˙ is the estimated elastic torque affected by the motor-side dynamics and it is calculated with nominal model parameters. In practice, the model differences will deteriorate the control performance. And normally it will be considered as a system disturbance which can be handled with disturbance observer method or parameter adaptive method.


## 2.1. Link-Side Dynamic of Robot

The variable stiffness actuated robot can be treated as complaint joints with stiffness variation. The general dynamic model of rigid robot with compliant joints and dynamic payload is recalled here. Without loss of generality, consider the variable stiffness actuated robot as an open kinematic chain havingN+1 rigid links, interconnected by N joints, and actuated by M electrical motors. The link-side coordinate q∈RN describes the joint angles and the motor-side coordinate θ∈RM denotes the electric motor angles. Obviously, the variable stiffness actuated robot is redundantly driven by electrical motors. These motors can be named as position motor and stiffness motor which are in charge of the position modulation and for output stiffness modulation of robot link, respectively. Therefore, the motor-side coordinate is defined as θ=[θpos,θstiff]. And θpos∈RN are the motor angles for position motor and θstiff∈RM-N are the motor angles of stiffness motor.Then, the kinetic energy, gravitational potential energy, and elastic potential energy are defined using the same formalism in [25]. Following the Lagrangian approach, the equation of motion is derived and represented as follows.( 1 ) Kinetic energy: the total kinetic energy of the rigid body with dynamic payload is given by(1)Ek=∑i=1NKi+12mLv0L22,where dynamic payload is denoted asmL and is considered as lumped mass which is fixed at the end-effector of robot. v0L is the velocity vector of the payload which is relative to the base coordinate. Ki represents the kinetic energy of each rigid link.( 2 ) Gravitational potential energy: The gravitational potential energy depends on the configuration of the rigid body dynamics. And the center of the mass of the actuators is assumed that will not change due to the self-rotation of the rotors and the gears.(2)UG=∑i=0N-mig·r0ci-mLg·r0cL,wheremi is the mass of robot’s i-th rigid link and r0ci is the mass center vector of i-th rigid link relative to the base coordinate. r0cL is the mass center vector relative to the base coordinate of dynamic payload which is defined above.( 3 ) Elastic potential energy: Due to the elastic element in variable stiffness actuators, the conservative potential energy can be represented as follows:(3)UE=∑i=0N∫τiθstiff,φθ,qdφ,whereτ(θstiff,φ)∈RN=r(θstiff)F(φ,θstiff) represents the general form of the elastic torque vector. r∈RN×p denotes the position matrix of the moment arm and F∈Rp denotes the corresponding forces due to the elastic elements (normally is the spring). Based on Spong’s model of flexible joint of robot, the variable stiffness actuator could be considered as complaint actuator with virtual torsional spring and then the deformation of spring is represented as φ(θ,q) and simplified as φ in the rest of the paper. In general, there are two different ways to modulate the output stiffness of the variable stiffness actuator: one way is through the adjustment of the spring preload [2] and the other way is through modulating the transmission ratio between input and output torque, such as lever principle [26].Following the Lagrangian approach, the equation of motion is derived and represented as follows:(4)L=EK-UG-UE(5)ddt∂L∂q˙-∂L∂q=0.Then, the dynamic model of variable stiffness actuated robot with dynamic payload is represented as(6)Mq,mLq¨+Cq,q˙,mLq˙+Dqq˙+Gq,mL=τEφ,θstiff,whereM(q,mL)∈RN×N is the symmetric and positive definite inertia matrix of the rigid body dynamics. C(q,q˙,mL)∈RN×N is the matrix that represents the Coriolis term, Dqq˙∈RN are the forces due to viscous friction, G(q,mL)∈RN are the gravitational forces and influenced only by the configuration of the robot, and τE(φ,θstiff)∈RN are the elastic joint torques that affect the rigid body dynamics.


## 2.2. Motor-Side Dynamic Model of Robot

The redundant motors are responsible for the output position and stiffness of robot links. Through the control of motor angles, we can obtain the stiffness and position of robot links simultaneously. Considering the motor-side dynamic model which is coupled with link-side dynamics will complicate the solution of the optimal control problem. Fortunately, for the servomotors with a high transmission ratio or high-gain feedback position controllers [27], the motor-side dynamics can be considered as decoupled from the link side and are expressed as follows:(7)JMθ¨+Dθθ˙=τEφ,θstiff+τM,whereJM∈RM×M represents the inertia of the motors which is a constant and diagonal matrix. Dθθ˙∈RM represents the viscous friction forces in the drive train. τM∈RM is the output torque of servomotors. Then, the classical second-order position controller of servomotors can be designed as(8)τM=τ^q,θ-ξθ˙-κθ-θd,whereθd∈RM is the reference position of motors for the control loop. ξ=diag⁡{2αi}∈RM×M, κ=diag⁡{αi2}∈RM×M are user defined gains as critically damped controller which leads to a fast response without overshooting control effort. And τ^=-τ^E(φ,θstiff)+D^θθ˙ is the estimated elastic torque affected by the motor-side dynamics and it is calculated with nominal model parameters. In practice, the model differences will deteriorate the control performance. And normally it will be considered as a system disturbance which can be handled with disturbance observer method or parameter adaptive method.


## 3. Optimal Control Problem Formulation

In this section, the weight lifting problem of variable stiffness actuated robot is formulated as a constrained optimal control problem in detail.


### 3.1. Optimal Objective Function

To maximum the payload of variable stiffness actuated robot, the optimal objective function is expressed as(9)J=-ωLmL+ωtimeT+∫0TτETωenergyτEdt,whereT is defined as duration of task execution. ωL∈R is the user-defined constant which is weight factor of dynamic payload. ωtime∈R is weight factor of execution time. ωenergy∈RN×N is user-defined constant matrix which is weight factor of the energy cost. The integral of the square of the elastic torque which represents the magnitude of the energy consumption [20, 28, 29] is utilized. We consider that the optimal control policy of the weight lifting problem is to obtain payload maximization along a point-to-point motion in a minimum time and energy cost under dynamic and kinematic constraints. Therefore, through choosing different weight factors among execution time, energy cost, and maximum payload, we can exploit the optimal policy like human-like performance in weight lifting games.


### 3.2. System Dynamics

The state-space representation of variable stiffness actuated robot is presented as follows. First of all, the state variable is defined as(10)x=x1,x2,x3,x4T=q,q˙,σ,σ˙T,where x∈R2(N+M). And the state variable σ is defined as σ=[φ,θstiff]. The initial state is defined as x0 and the target state is defined as xf.(11)x0=x0,xT=xf.And the control variable is defined as(12)u=u1,u2=θstiffd,φd,where u∈RM. And φd,θstiffd represent reference deformation and motor angles of the control loop at the motor-side dynamic, respectively. Then, the system dynamics can be represented as(13)x˙=fx,u,t=x2Mx1,mL-1τEx1,x3-Cx1,x2,mL-Dx2-Gx1,mLx4-ξx4-κx3+κu.


### 3.3. State Constraints

Due to the physical limits on the deformation of the elastic elements, like the spring deformation, nonlinear state inequality constraints of variable stiffness actuators should be satisfied as follows:(14)φmin≤φq,θ≤φmax,whereφ(q,θ)∈RM is defined as deformation of virtual torsional spring related to link angles and motor positions. And φmin and φmax represent the lower and upper bounds of deformation of the virtual spring, respectively. However, particularly, in the case of variable stiffness actuated robot, due to the physical restriction of the mechanical structure, the violation of the deformation constraints in (14) would cause permanent damage to the structure. Therefore, the deformation constraints should be strictly satisfied for the real-world application in advance.What is more, the link trajectory is subject to the inequality constraints posed by the range and speed limit at link side. Therefore, the admissible set of states subject to these constraints can be defined by(15)qmin≤q≤qmaxq˙min≤q˙≤q˙max,whereq˙min and q˙max are the lower and upper bounds of link angular velocity. qmin and qmax are the lower and upper bounds of link angles.Based on the physical limit of servomotors and the mechanical structure, the motor trajectories are subject to inequality constraints due to the range, rate, and acceleration limits. The admissible set of trajectories subject to these constraints can be defined by(16)θmin≤θ≤θmaxθ˙min≤θ˙≤θ˙maxθ¨min≤θ¨≤θ¨max,whereθ˙min and θ˙max are the lower and upper bounds of motors angular velocity. θ¨min and θ¨max are the lower and upper bounds of motor acceleration.In summary, the state constraints can be collected as an admissible set as follows:(17)H=x∈R2M+N:φmin≤φq,θ≤φmaxqmin≤q≤qmaxq˙min≤q˙≤q˙maxθmin≤θ≤θmaxθ˙min≤θ˙≤θ˙maxθ¨min≤θ¨≤θ¨max.


### 3.4. Optimal Control Problem

Therefore, the weight lifting problem with variable stiffness actuated robot is expressed as a constrained optimal control problem as follows:(18)minT,u,mLJ=-ωLmL+ωtimeT+∫0TτETωenergyτEdts.t.x˙=fx,u,t,x0=x0,xT=xf,x∈H,u∈U.The optimal control problem above can be solved with different types of optimization method, like direct method or indirect method. In this paper, a subset of direct collocation methods called Gauss pseudospectral method is employed.


## 3.1. Optimal Objective Function

To maximum the payload of variable stiffness actuated robot, the optimal objective function is expressed as(9)J=-ωLmL+ωtimeT+∫0TτETωenergyτEdt,whereT is defined as duration of task execution. ωL∈R is the user-defined constant which is weight factor of dynamic payload. ωtime∈R is weight factor of execution time. ωenergy∈RN×N is user-defined constant matrix which is weight factor of the energy cost. The integral of the square of the elastic torque which represents the magnitude of the energy consumption [20, 28, 29] is utilized. We consider that the optimal control policy of the weight lifting problem is to obtain payload maximization along a point-to-point motion in a minimum time and energy cost under dynamic and kinematic constraints. Therefore, through choosing different weight factors among execution time, energy cost, and maximum payload, we can exploit the optimal policy like human-like performance in weight lifting games.


## 3.2. System Dynamics

The state-space representation of variable stiffness actuated robot is presented as follows. First of all, the state variable is defined as(10)x=x1,x2,x3,x4T=q,q˙,σ,σ˙T,where x∈R2(N+M). And the state variable σ is defined as σ=[φ,θstiff]. The initial state is defined as x0 and the target state is defined as xf.(11)x0=x0,xT=xf.And the control variable is defined as(12)u=u1,u2=θstiffd,φd,where u∈RM. And φd,θstiffd represent reference deformation and motor angles of the control loop at the motor-side dynamic, respectively. Then, the system dynamics can be represented as(13)x˙=fx,u,t=x2Mx1,mL-1τEx1,x3-Cx1,x2,mL-Dx2-Gx1,mLx4-ξx4-κx3+κu.


## 3.3. State Constraints

Due to the physical limits on the deformation of the elastic elements, like the spring deformation, nonlinear state inequality constraints of variable stiffness actuators should be satisfied as follows:(14)φmin≤φq,θ≤φmax,whereφ(q,θ)∈RM is defined as deformation of virtual torsional spring related to link angles and motor positions. And φmin and φmax represent the lower and upper bounds of deformation of the virtual spring, respectively. However, particularly, in the case of variable stiffness actuated robot, due to the physical restriction of the mechanical structure, the violation of the deformation constraints in (14) would cause permanent damage to the structure. Therefore, the deformation constraints should be strictly satisfied for the real-world application in advance.What is more, the link trajectory is subject to the inequality constraints posed by the range and speed limit at link side. Therefore, the admissible set of states subject to these constraints can be defined by(15)qmin≤q≤qmaxq˙min≤q˙≤q˙max,whereq˙min and q˙max are the lower and upper bounds of link angular velocity. qmin and qmax are the lower and upper bounds of link angles.Based on the physical limit of servomotors and the mechanical structure, the motor trajectories are subject to inequality constraints due to the range, rate, and acceleration limits. The admissible set of trajectories subject to these constraints can be defined by(16)θmin≤θ≤θmaxθ˙min≤θ˙≤θ˙maxθ¨min≤θ¨≤θ¨max,whereθ˙min and θ˙max are the lower and upper bounds of motors angular velocity. θ¨min and θ¨max are the lower and upper bounds of motor acceleration.In summary, the state constraints can be collected as an admissible set as follows:(17)H=x∈R2M+N:φmin≤φq,θ≤φmaxqmin≤q≤qmaxq˙min≤q˙≤q˙maxθmin≤θ≤θmaxθ˙min≤θ˙≤θ˙maxθ¨min≤θ¨≤θ¨max.


## 3.4. Optimal Control Problem

Therefore, the weight lifting problem with variable stiffness actuated robot is expressed as a constrained optimal control problem as follows:(18)minT,u,mLJ=-ωLmL+ωtimeT+∫0TτETωenergyτEdts.t.x˙=fx,u,t,x0=x0,xT=xf,x∈H,u∈U.The optimal control problem above can be solved with different types of optimization method, like direct method or indirect method. In this paper, a subset of direct collocation methods called Gauss pseudospectral method is employed.


## 4. Numeric Solution with Gauss Pseudospectral Method

In this section, the Gauss pseudospectral method is employed to solve the optimal control problem formulated for the weight lifting problem with variable stiffness actuated robot. This method utilizes global polynomial approximations at a set of Gauss collocation points for the control and state variables. Thus, the continuous optimal control problem can be transformed into a discrete nonlinear programming problem (NLP), which can be solved with plenty of well-developed solvers. Although there are different types of pseudospectral method featured with different set of collocation points, like Legendre-Gauss (LG), Legendre-Gauss-Radau (LGR), and Legendre-Gauss-Lobatto (LGL) points, the set of Legendre-Gauss points is utilized in this paper for the accurate approximation at start and end nodes [21].


### 4.1. Time Domain Transformation

Because the Gauss collocation points are interior to the interval[-1,1], the original time interval t∈[t0,tf] has to be scaled to the time interval s∈[-1,1] by(19)s=2t-tf-t0tf-t0.Based on Gauss’s theory, the polynomial approximation with the Guass collocation points can maintain approximation accuracy at the interpolation nodes. TheN-th Legendre-Gauss points [s1,s2,…,sN] and s0=-1 are chosen and defined as interpolation nodes [s0,s1,s2,…,sN]. (20)-1=s0<s1<s2<s3<⋯<sN<1,where the Legendre-Gauss points[s1,s2,…,sN] are the root of N-th Legendre polynomials derived from PN(s)=0, namely,(21)PNs=12NN!dNdsNs2-1N=0.Based on the stated above, the discretization of system dynamics, state, and control variables is presented as follows.


### 4.2. The State and Control Approximation

To approximate the state and control variables with the set of Legendre-Gauss points above, a basis function of(N+1)-th Lagrange interpolating polynomials and the state variables are approximated as follows:(22)xs≈Xs=∑i=0NxsiLis,where the Lagrange interpolation functions are(23)Lis=∏j=0,j≠iNs-sjsi-sj.Meanwhile, the control variables are also approximated usingN-th Lagrange interpolating polynomials because the number of the interpolation nodes is N. Then, the approximation for control variables can be expressed as follows:(24)us≈Us=∑i=0NusiLi∗s,where the Lagrange interpolation functions are(25)Li∗s=∏j=1,j≠iNs-sjsi-sj.


### 4.3. The System Dynamic Approximation

Differentiating system dynamic, we obtain(26)x˙s≈X˙s=∑i=0NxsiL˙is.In the Guass pseudospectral method, the dynamics are collocated at the N Legendre–Guass point. The derivative of each Lagrange polynomial at the LG points can be represented in a differential approximation matrixDki. The differential matrix can be represented as follows:(27)
D
k
i
=
L
˙
i
s
k
=
d
L
i
d
s
s
k
=
1
+
s
k
P
˙
k
s
k
+
P
k
s
k
s
k
-
s
i
1
+
s
i
P
˙
k
s
k
+
P
k
s
k
,
i
≠
k
1
+
s
k
P
¨
k
s
k
+
2
P
˙
k
s
k
2
1
+
s
i
P
˙
k
s
k
+
P
k
s
k
,
i
=
k
,wherek=1,…,N and i=0,…,N. Thus, the dynamic differential equation constraints are transformed into discrete constraints:(28)∑i=0NDkiXsi-tf-t02fXsk,Usk,sk=0.


### 4.4. The Terminal State Approximation

However, the terminal node is not included in the approximate expression of state variables. Then, the terminal state is obtained as follows:(29)xtf=x0+∫-11fxs,us,sds.


### 4.5. Optimal Objective Function Approximation

When the integral item of objective function is approximated via the Gauss integral, the GPM approximation of the objective function can be obtained as(30)J=ϕX0,t0,Xf,tf+tf-t02∑i=1NωifXi,Ui,si,whereϕ(X0,t0,Xf,tf)∈R is the terminal performance function. ωi are the Gauss weights which are calculated as follows:(31)ωi=21-si2Pisi2.To sum up, the continuous control and state variables are discretized at the Guass collocation points. And then the system dynamics and optimal objective function are transformed and represented with the discretized variables. Therefore, for the optimal control problem, we redefine the new decision variables as(32)z=Xi,Ui,t0,tf,mL.Thus, the original optimal control problem of variable stiffness actuated robot can be transformed as a general nonlinear programming problem as follows:(33)minzJzs.t.gz=0,hz≤0,whereg(z),h(z) denote the equality and inequality constraints, respectively. The converted nonlinear programming problem is a high-dimensional and sparse problem. For the solution of the general nonlinear programming problem above, there are plenty of choices of mature numerical sparse solvers, like IPOPT in this paper.


## 4.1. Time Domain Transformation

Because the Gauss collocation points are interior to the interval[-1,1], the original time interval t∈[t0,tf] has to be scaled to the time interval s∈[-1,1] by(19)s=2t-tf-t0tf-t0.Based on Gauss’s theory, the polynomial approximation with the Guass collocation points can maintain approximation accuracy at the interpolation nodes. TheN-th Legendre-Gauss points [s1,s2,…,sN] and s0=-1 are chosen and defined as interpolation nodes [s0,s1,s2,…,sN]. (20)-1=s0<s1<s2<s3<⋯<sN<1,where the Legendre-Gauss points[s1,s2,…,sN] are the root of N-th Legendre polynomials derived from PN(s)=0, namely,(21)PNs=12NN!dNdsNs2-1N=0.Based on the stated above, the discretization of system dynamics, state, and control variables is presented as follows.


## 4.2. The State and Control Approximation

To approximate the state and control variables with the set of Legendre-Gauss points above, a basis function of(N+1)-th Lagrange interpolating polynomials and the state variables are approximated as follows:(22)xs≈Xs=∑i=0NxsiLis,where the Lagrange interpolation functions are(23)Lis=∏j=0,j≠iNs-sjsi-sj.Meanwhile, the control variables are also approximated usingN-th Lagrange interpolating polynomials because the number of the interpolation nodes is N. Then, the approximation for control variables can be expressed as follows:(24)us≈Us=∑i=0NusiLi∗s,where the Lagrange interpolation functions are(25)Li∗s=∏j=1,j≠iNs-sjsi-sj.


## 4.3. The System Dynamic Approximation

Differentiating system dynamic, we obtain(26)x˙s≈X˙s=∑i=0NxsiL˙is.In the Guass pseudospectral method, the dynamics are collocated at the N Legendre–Guass point. The derivative of each Lagrange polynomial at the LG points can be represented in a differential approximation matrixDki. The differential matrix can be represented as follows:(27)
D
k
i
=
L
˙
i
s
k
=
d
L
i
d
s
s
k
=
1
+
s
k
P
˙
k
s
k
+
P
k
s
k
s
k
-
s
i
1
+
s
i
P
˙
k
s
k
+
P
k
s
k
,
i
≠
k
1
+
s
k
P
¨
k
s
k
+
2
P
˙
k
s
k
2
1
+
s
i
P
˙
k
s
k
+
P
k
s
k
,
i
=
k
,wherek=1,…,N and i=0,…,N. Thus, the dynamic differential equation constraints are transformed into discrete constraints:(28)∑i=0NDkiXsi-tf-t02fXsk,Usk,sk=0.


## 4.4. The Terminal State Approximation

However, the terminal node is not included in the approximate expression of state variables. Then, the terminal state is obtained as follows:(29)xtf=x0+∫-11fxs,us,sds.


## 4.5. Optimal Objective Function Approximation

When the integral item of objective function is approximated via the Gauss integral, the GPM approximation of the objective function can be obtained as(30)J=ϕX0,t0,Xf,tf+tf-t02∑i=1NωifXi,Ui,si,whereϕ(X0,t0,Xf,tf)∈R is the terminal performance function. ωi are the Gauss weights which are calculated as follows:(31)ωi=21-si2Pisi2.To sum up, the continuous control and state variables are discretized at the Guass collocation points. And then the system dynamics and optimal objective function are transformed and represented with the discretized variables. Therefore, for the optimal control problem, we redefine the new decision variables as(32)z=Xi,Ui,t0,tf,mL.Thus, the original optimal control problem of variable stiffness actuated robot can be transformed as a general nonlinear programming problem as follows:(33)minzJzs.t.gz=0,hz≤0,whereg(z),h(z) denote the equality and inequality constraints, respectively. The converted nonlinear programming problem is a high-dimensional and sparse problem. For the solution of the general nonlinear programming problem above, there are plenty of choices of mature numerical sparse solvers, like IPOPT in this paper.


## 5. Case Study of 2-DOF Variable Stiffness Actuated Robot

In this section, a model of 2-DOF (Degree Of Freedom) variable stiffness actuated robot is built to verify the proposed method for the weight lifting problem and its schematic drawing is depicted in Figure1. A variable stiffness actuator based on roller-cam mechanism is utilized to drive the link of the robot in Figure 2. The design and stress analysis of this actuator have been presented in detail in the authors’ early work [19]. Nevertheless, the model of the actuator and the deformation constraints are recalled in this paper. Then, to show the effectiveness of the proposed method, the optimization results are compared with the results of iLQR approach. As shown in (30), the optimization results are highly relied on the choice of the weight factor in the objective function. Therefore, a three-dimensional parameter grid matrix has come up and the optimizations with different setups of objective function parameters on the grid matrix are conducted. Finally, how the variable stiffness actuated robot can perform better in the weight lifting problem is analyzed.Figure 1
Schematic drawing of two-link variable stiffness actuated robot and the physical parameters of dynamics are reported in Table1.Table 1
Simulation parameters.
Symbol Parameter Value/Unit ω t i m e weight factor of duration 0.5 ω e n e r g y weight factor of energy cost 1 e - 4 0 0 1 e - 4 ω L weight factor of payload 2 φ m a x maximum deformation of elastic element 0.2617 r a d κ stiffness coefficient in motor position controller 250 0 0 0 0 250 0 0 0 0 100 0 0 0 0 100   N / m ξ damping coefficient in motor position controller 50 0 0 0 0 50 0 0 0 0 10 0 0 0 0 10   N / ( m / s ) ρ c h a r characteristic coefficient of cam disk curve 12 ε the optimization tolerance 0.5 R the moment arm of elastic torque 33mm K spring stiffness 8000N/m τ m a x bounds of the maximum output torque [ - 5,5 ] N ⋅ m [ θ m i n s t i f f , θ m a x s t i f f ] bounds of stiffness motor angle [ 0,0.174 ] r a d  [ θ ˙ m i n s t i f f , θ ˙ m a x s t i f f ] bounds of stiffness motor velocity [ - 1.96,1.96 ] r a d / s  [ θ ¨ m i n s t i f f , θ ¨ m a x s t i f f ] bounds of stiffness motor acceleration [ - 10,10 ] r a d / s [ B 1 , B 2 ] viscous friction coefficient [ 0.01,0.01 ] N m s / r a d [ m 1 , m 2 ] mass of link1 and link2 [ 0.5,0.5 ] k g [ I 1 , I 2 ] inertia of link1 and link2 [ 0.5,0.5 ] k g ⋅ m 2 [ L 1 , L 2 ] length of link1 and link2 [ 0.5,0.5 ] m [ L c 1 , L c 2 ] the center of mass of link1 and link2 [ 0.5,0.5 ] m x 0 initial state [ - 1.57,0 , 0,0 , 0,0 , 0,0 , 0,0 , 0,0 ]  x f terminal state [ 1.57,0 , 0,0 , 0,0 , 0,0 , 0,0 , 0,0 ]Figure 2
Mechanical design of variable stiffness actuator based on roller-cam mechanism. Three pairs of rollers are directly connected to the output link and driven by position motor through harmonic reducer. The cam disks are axially pushed together by the linear spring. Then, the stiffness setup can be adjusted by the relatively rotation of the two cam disks through the stiffness motors.


### 5.1. Model of 2-DOF Variable Stiffness Actuated Robot

#### 5.1.1. System Dynamics

The vector of link-side angles describing the rotation of the two links is defined asq=[q1,q2] and the inertia matrix M(q)=m11m12m21m22, the Coriolis and normal inertial terms C(q,q˙)=[c1,c2], the viscous friction D(q˙)=[d1,d2], and the gravity terms G(q)=[g1,g2] are defined as follows:(34)m11=I1+m1Lc1+I2+m2L12+2L1Lc2cos⁡q2+Lc22+mLL12+2L1L2cos⁡q2+L22m12=m21=I2+m2Lc22+2L1Lc2cos⁡q2+mLL22+2L1L2cos⁡q2m22=I2+m2Lc22+mLL22c1=-m2Lc2+mLL2L1sin⁡q22q˙1q˙2+q˙22c2=m2Lc2+mLL2L1sin⁡q2q˙12d1=B1q˙1d2=B2q˙2g1=gm1Lc1+m2L-1+mLL1cos⁡q1+gm2Lc2+mLL2cos⁡q1+q2g2=gm2Lc2+mLL2cos⁡q1+q2,wheremi,Li,Lci,Ii,Bi,i=1,2 are, respectively, the mass, the length of the i-th link, the center of mass, and the moment of inertia of the same link.


#### 5.1.2. Actuation Constraints

In case of the robot illustrated in Figure1, the deformation of virtual torsional spring is defined as φ∈{R2:φi=qi-θi+2,i=1,2}. θstiff=[θ2,θ4]∈R2 represents the stiffness setup angles. From the stress analysis in the early work [19] of the variable stiffness actuator in Figure 2, the output elastic torque is formulated as follows:(35)τE1=τmaxeρchar-φmax+θstiffeρcharφ-e-ρcharφ,whereτE1=[τ1,τ2]∈R2. τmax∈R2 is a user-defined maximum output torque vector when the maximum deformation of spring is achieved. ρchar is characteristic constant of the designed cam disk profile. Based on the mechanism of the variable stiffness actuator, the deformation of the virtual spring is physically limited by the compression length of linear spring. And the maximum of the deformation is also related to the stiffness setup. Then, the relation between deformation constraints and stiffness setup is shown as follows:(36)-φmax+θstiff≤φ≤φmax-θstiff.As stated in Section4, the actuator constraints are discretized and treated as general inequality constraints. However, the efficiency of the numerical solvers is very sensitive to the number of the constraints. Therefore, a reduction in the number and complexity of the constraints is often preferred in practice. Then, the inequality constraints are embedded into the system dynamics through the state space transformation, instead of being treated as penalty terms in the objective function [19].The transformation of the stateσ=[φ,θstiff] and the new state variables σ¯ are stated as follows:(37)Δ=σ¯∈0,1:σ¯i=σi-σminiσmaxi-σmini,i=1,…,M.Then, substitute the new state variablesσ[i]=σmin[i]+σmax[i]-σmin[i]  σ¯[i] into the system dynamics (13). The state transformation makes the setting of the initial guess easier and therefore improves the success rate of the optimization dramatically.


### 5.2. Optimization Results

In this section, to verify the proposed method, the optimization simulations of the weight lifting problem along a predefined point-to-point motion on variable stiffness actuated robot are carried out. As shown in (30), the maximum payload of robot is treated as negative term in the objective function. To verify the effectiveness of the payload maximization term, optimizations with and without payload term are conducted with the proposed method based on Guass pseudospectral method. After payload optimization in this paper, the maximum payload of the robot is improved to 3.8 kg. And the trajectory of the robot link and stiffness motor angles are depicted in Figures 5 and 6, respectively. As illustrated in Figure 6, it is noteworthy that, to maximize the payload which robot can lift out, the optimal control policy takes advantage of the gravity. Through the swing of robot links, it costs less input energy to lift the maximum weight while the lifting energy is from the gravity potential energy rather than directly from the output torque of the actuators. What is more, the optimal results show the output torques have been saturated more sufficiently in Figure 7. Furthermore, to confirm the effectiveness of the proposed method, the same optimal control problem is also solved with the iLQR approach. However, the payload parameter of the system dynamics cannot be approximated directly under the iLQR framework; then, a two-dimensional grid of duration and payload has come up. Thus, we can only solve the energy optimal control problem repeatedly with fixed duration and payload at the predefined parameters grid and the solution with maximum payload and minimum time will be the optimal one. The optimization results with iLQR approach are depicted in Figures 5, 3, 6, and 7. The optimization trajectory with proposed method is much more smooth than the results with iLQR approach. It is because the control and state variables have been approximated with global polynomial interpolation function with Guass pseudospectral method. And under iLQR framework, only the suboptimal solution will be found instead of the global ones. Therefore, the smooth trajectory obtained with the proposed method will be more easily treated as reference signals of the physical servomotor system.Figure 3
Optimal deformation trajectories of variable stiffness actuated robot with GPM and iLQR. The blue dotted lines and the red dashed lines are the lower and upper bounds defined in (36) while they are not violated during the entire point-to-point motion. For comparison in the same time scale, the path parameter {s∈R,0≤s≤1} is introduced to scale different execution times into the nominal parameter space. And in the rest of paper, the comparisons among different actuators are also transformed to the parameter space.
(a)
Actuator deformation with Gauss pseudospectral method (GPM) (b)
Actuator deformation with iterative linear quadratic regulator (iLQR)The comparison of the optimization results have been shown in Table2. The target accuracy of iLQR approach is relatively worse than the proposed method. However, the computation time of the optimization is faster than the proposed method due to the suboptimal solution at each optimization iteration. Therefore, the iLQR approach has more potential to solve the optimal control problem online. In practice, the iLQR approach is more sensitive to the initial guess than the proposed method. Therefore, to set a suitable initial guess for the optimization process, the random initial guess has to be utilized. As shown in Table 2, the larger maximum payload is achieved with the proposed method while the motion duration is slightly smaller than the results obtained by iLQR method. Based on the optimization results, the proposed method is more effective than the iLQR method for offline trajectory optimization.Table 2
Comparison of optimization results.
Maximum payload(kg) Motion duration(s) Computation time(s) Target accuracy1 G P M 2 0 1.72 95.068 1e-6 iLQR 2.5 5 42.9 0.1647 G P M h s 3 3.9 4.91 151.2166 1 e - 6 G P M l s 4 3.8 4.95 115.842 1 e - 6 1 Target accuracy computed as xf∗-xf2 which reflects the deviation from the terminal state 2  ωL is set to zero 3 The stiffness motor ranges are set as higher stiffness which is [0.087,0.174]rad. 4 The ranges of stiffness motors are set as lower stiffness which is [0,0.087]rad.Nevertheless, as shown in (30), the optimization results highly relied on the choice of the weight factors in the objective function. To study the trade-off of the time, energy cost, and maximum payload in the weight lifting problem with robot, a three-dimensional parameter grid has come up and the optimizations with different setups of weight factors at the grid are conducted. The results are illustrated in Figure 4. The weight factors of time, energy cost, and payload ranged from 0 to 1, 1 and 2, respectively. Then, the three-dimensional grid matrix is meshed with fixed interval. The optimization results show that it is not necessary to choose the largest weight factors to get maximum payload in a short time. However, the sweet area of better desired performance is still around the large weight factors setups.Figure 4
The optimization results with different setups of weight factors.
(a)
The results of maximum payload at the weight factors grid (b)
The results of maximum time at the weight factors gridFigure 5
The motor trajectories of stiffness modulation.
(a)
Motor anglesθ2 for the stiffness modulation (b)
Motor anglesθ4 for the stiffness modulationFigure 6
The link-side trajectory.
(a)
Joint angleq1 (b)
Joint angleq2Figure 7
The joint torques.
(a)
Joint torque of link1 (b)
Joint torque of link2In summary, the effectiveness of the proposed method has been proved. However, how the variable stiffness actuated robot can perform better in the weight lifting problem still should be answered. With lower and higher setup of stiffness of the variable stiffness actuators, simulations are carried out to find the influence of the stiffness variation on weight lifting motion. As shown in Table2, there is no big influence on the maximum payload and minimum time with high or lower stiffness setup. However, the torque limits of each actuator have direct influence on the maximum payload. Higher torque limit will lead to more output energy for the weight lifting. On the other side, at the design of variable stiffness actuators for the weight lifting problem, extending the torque saturation is preferred than the ability of stiffness variation with variable stiffness actuated robot. The guideline is similar with the athletes’ training program in which the pursuit of more muscle strength is more important than the skills for the weight lifting sports. And it is interesting to notice that the property of the variable stiffness actuators has to meet the core requirement in the different applications, like stiffness variation for the explosive movement or periodic movement and larger torque saturation in the same size and input power for the weight lifting movement.


## 5.1. Model of 2-DOF Variable Stiffness Actuated Robot

### 5.1.1. System Dynamics

The vector of link-side angles describing the rotation of the two links is defined asq=[q1,q2] and the inertia matrix M(q)=m11m12m21m22, the Coriolis and normal inertial terms C(q,q˙)=[c1,c2], the viscous friction D(q˙)=[d1,d2], and the gravity terms G(q)=[g1,g2] are defined as follows:(34)m11=I1+m1Lc1+I2+m2L12+2L1Lc2cos⁡q2+Lc22+mLL12+2L1L2cos⁡q2+L22m12=m21=I2+m2Lc22+2L1Lc2cos⁡q2+mLL22+2L1L2cos⁡q2m22=I2+m2Lc22+mLL22c1=-m2Lc2+mLL2L1sin⁡q22q˙1q˙2+q˙22c2=m2Lc2+mLL2L1sin⁡q2q˙12d1=B1q˙1d2=B2q˙2g1=gm1Lc1+m2L-1+mLL1cos⁡q1+gm2Lc2+mLL2cos⁡q1+q2g2=gm2Lc2+mLL2cos⁡q1+q2,wheremi,Li,Lci,Ii,Bi,i=1,2 are, respectively, the mass, the length of the i-th link, the center of mass, and the moment of inertia of the same link.


### 5.1.2. Actuation Constraints

In case of the robot illustrated in Figure1, the deformation of virtual torsional spring is defined as φ∈{R2:φi=qi-θi+2,i=1,2}. θstiff=[θ2,θ4]∈R2 represents the stiffness setup angles. From the stress analysis in the early work [19] of the variable stiffness actuator in Figure 2, the output elastic torque is formulated as follows:(35)τE1=τmaxeρchar-φmax+θstiffeρcharφ-e-ρcharφ,whereτE1=[τ1,τ2]∈R2. τmax∈R2 is a user-defined maximum output torque vector when the maximum deformation of spring is achieved. ρchar is characteristic constant of the designed cam disk profile. Based on the mechanism of the variable stiffness actuator, the deformation of the virtual spring is physically limited by the compression length of linear spring. And the maximum of the deformation is also related to the stiffness setup. Then, the relation between deformation constraints and stiffness setup is shown as follows:(36)-φmax+θstiff≤φ≤φmax-θstiff.As stated in Section4, the actuator constraints are discretized and treated as general inequality constraints. However, the efficiency of the numerical solvers is very sensitive to the number of the constraints. Therefore, a reduction in the number and complexity of the constraints is often preferred in practice. Then, the inequality constraints are embedded into the system dynamics through the state space transformation, instead of being treated as penalty terms in the objective function [19].The transformation of the stateσ=[φ,θstiff] and the new state variables σ¯ are stated as follows:(37)Δ=σ¯∈0,1:σ¯i=σi-σminiσmaxi-σmini,i=1,…,M.Then, substitute the new state variablesσ[i]=σmin[i]+σmax[i]-σmin[i]  σ¯[i] into the system dynamics (13). The state transformation makes the setting of the initial guess easier and therefore improves the success rate of the optimization dramatically.


## 5.1.1. System Dynamics

The vector of link-side angles describing the rotation of the two links is defined asq=[q1,q2] and the inertia matrix M(q)=m11m12m21m22, the Coriolis and normal inertial terms C(q,q˙)=[c1,c2], the viscous friction D(q˙)=[d1,d2], and the gravity terms G(q)=[g1,g2] are defined as follows:(34)m11=I1+m1Lc1+I2+m2L12+2L1Lc2cos⁡q2+Lc22+mLL12+2L1L2cos⁡q2+L22m12=m21=I2+m2Lc22+2L1Lc2cos⁡q2+mLL22+2L1L2cos⁡q2m22=I2+m2Lc22+mLL22c1=-m2Lc2+mLL2L1sin⁡q22q˙1q˙2+q˙22c2=m2Lc2+mLL2L1sin⁡q2q˙12d1=B1q˙1d2=B2q˙2g1=gm1Lc1+m2L-1+mLL1cos⁡q1+gm2Lc2+mLL2cos⁡q1+q2g2=gm2Lc2+mLL2cos⁡q1+q2,wheremi,Li,Lci,Ii,Bi,i=1,2 are, respectively, the mass, the length of the i-th link, the center of mass, and the moment of inertia of the same link.


## 5.1.2. Actuation Constraints

In case of the robot illustrated in Figure1, the deformation of virtual torsional spring is defined as φ∈{R2:φi=qi-θi+2,i=1,2}. θstiff=[θ2,θ4]∈R2 represents the stiffness setup angles. From the stress analysis in the early work [19] of the variable stiffness actuator in Figure 2, the output elastic torque is formulated as follows:(35)τE1=τmaxeρchar-φmax+θstiffeρcharφ-e-ρcharφ,whereτE1=[τ1,τ2]∈R2. τmax∈R2 is a user-defined maximum output torque vector when the maximum deformation of spring is achieved. ρchar is characteristic constant of the designed cam disk profile. Based on the mechanism of the variable stiffness actuator, the deformation of the virtual spring is physically limited by the compression length of linear spring. And the maximum of the deformation is also related to the stiffness setup. Then, the relation between deformation constraints and stiffness setup is shown as follows:(36)-φmax+θstiff≤φ≤φmax-θstiff.As stated in Section4, the actuator constraints are discretized and treated as general inequality constraints. However, the efficiency of the numerical solvers is very sensitive to the number of the constraints. Therefore, a reduction in the number and complexity of the constraints is often preferred in practice. Then, the inequality constraints are embedded into the system dynamics through the state space transformation, instead of being treated as penalty terms in the objective function [19].The transformation of the stateσ=[φ,θstiff] and the new state variables σ¯ are stated as follows:(37)Δ=σ¯∈0,1:σ¯i=σi-σminiσmaxi-σmini,i=1,…,M.Then, substitute the new state variablesσ[i]=σmin[i]+σmax[i]-σmin[i]  σ¯[i] into the system dynamics (13). The state transformation makes the setting of the initial guess easier and therefore improves the success rate of the optimization dramatically.


## 5.2. Optimization Results

In this section, to verify the proposed method, the optimization simulations of the weight lifting problem along a predefined point-to-point motion on variable stiffness actuated robot are carried out. As shown in (30), the maximum payload of robot is treated as negative term in the objective function. To verify the effectiveness of the payload maximization term, optimizations with and without payload term are conducted with the proposed method based on Guass pseudospectral method. After payload optimization in this paper, the maximum payload of the robot is improved to 3.8 kg. And the trajectory of the robot link and stiffness motor angles are depicted in Figures 5 and 6, respectively. As illustrated in Figure 6, it is noteworthy that, to maximize the payload which robot can lift out, the optimal control policy takes advantage of the gravity. Through the swing of robot links, it costs less input energy to lift the maximum weight while the lifting energy is from the gravity potential energy rather than directly from the output torque of the actuators. What is more, the optimal results show the output torques have been saturated more sufficiently in Figure 7. Furthermore, to confirm the effectiveness of the proposed method, the same optimal control problem is also solved with the iLQR approach. However, the payload parameter of the system dynamics cannot be approximated directly under the iLQR framework; then, a two-dimensional grid of duration and payload has come up. Thus, we can only solve the energy optimal control problem repeatedly with fixed duration and payload at the predefined parameters grid and the solution with maximum payload and minimum time will be the optimal one. The optimization results with iLQR approach are depicted in Figures 5, 3, 6, and 7. The optimization trajectory with proposed method is much more smooth than the results with iLQR approach. It is because the control and state variables have been approximated with global polynomial interpolation function with Guass pseudospectral method. And under iLQR framework, only the suboptimal solution will be found instead of the global ones. Therefore, the smooth trajectory obtained with the proposed method will be more easily treated as reference signals of the physical servomotor system.Figure 3
Optimal deformation trajectories of variable stiffness actuated robot with GPM and iLQR. The blue dotted lines and the red dashed lines are the lower and upper bounds defined in (36) while they are not violated during the entire point-to-point motion. For comparison in the same time scale, the path parameter {s∈R,0≤s≤1} is introduced to scale different execution times into the nominal parameter space. And in the rest of paper, the comparisons among different actuators are also transformed to the parameter space.
(a)
Actuator deformation with Gauss pseudospectral method (GPM) (b)
Actuator deformation with iterative linear quadratic regulator (iLQR)The comparison of the optimization results have been shown in Table2. The target accuracy of iLQR approach is relatively worse than the proposed method. However, the computation time of the optimization is faster than the proposed method due to the suboptimal solution at each optimization iteration. Therefore, the iLQR approach has more potential to solve the optimal control problem online. In practice, the iLQR approach is more sensitive to the initial guess than the proposed method. Therefore, to set a suitable initial guess for the optimization process, the random initial guess has to be utilized. As shown in Table 2, the larger maximum payload is achieved with the proposed method while the motion duration is slightly smaller than the results obtained by iLQR method. Based on the optimization results, the proposed method is more effective than the iLQR method for offline trajectory optimization.Table 2
Comparison of optimization results.
Maximum payload(kg) Motion duration(s) Computation time(s) Target accuracy1 G P M 2 0 1.72 95.068 1e-6 iLQR 2.5 5 42.9 0.1647 G P M h s 3 3.9 4.91 151.2166 1 e - 6 G P M l s 4 3.8 4.95 115.842 1 e - 6 1 Target accuracy computed as xf∗-xf2 which reflects the deviation from the terminal state 2  ωL is set to zero 3 The stiffness motor ranges are set as higher stiffness which is [0.087,0.174]rad. 4 The ranges of stiffness motors are set as lower stiffness which is [0,0.087]rad.Nevertheless, as shown in (30), the optimization results highly relied on the choice of the weight factors in the objective function. To study the trade-off of the time, energy cost, and maximum payload in the weight lifting problem with robot, a three-dimensional parameter grid has come up and the optimizations with different setups of weight factors at the grid are conducted. The results are illustrated in Figure 4. The weight factors of time, energy cost, and payload ranged from 0 to 1, 1 and 2, respectively. Then, the three-dimensional grid matrix is meshed with fixed interval. The optimization results show that it is not necessary to choose the largest weight factors to get maximum payload in a short time. However, the sweet area of better desired performance is still around the large weight factors setups.Figure 4
The optimization results with different setups of weight factors.
(a)
The results of maximum payload at the weight factors grid (b)
The results of maximum time at the weight factors gridFigure 5
The motor trajectories of stiffness modulation.
(a)
Motor anglesθ2 for the stiffness modulation (b)
Motor anglesθ4 for the stiffness modulationFigure 6
The link-side trajectory.
(a)
Joint angleq1 (b)
Joint angleq2Figure 7
The joint torques.
(a)
Joint torque of link1 (b)
Joint torque of link2In summary, the effectiveness of the proposed method has been proved. However, how the variable stiffness actuated robot can perform better in the weight lifting problem still should be answered. With lower and higher setup of stiffness of the variable stiffness actuators, simulations are carried out to find the influence of the stiffness variation on weight lifting motion. As shown in Table2, there is no big influence on the maximum payload and minimum time with high or lower stiffness setup. However, the torque limits of each actuator have direct influence on the maximum payload. Higher torque limit will lead to more output energy for the weight lifting. On the other side, at the design of variable stiffness actuators for the weight lifting problem, extending the torque saturation is preferred than the ability of stiffness variation with variable stiffness actuated robot. The guideline is similar with the athletes’ training program in which the pursuit of more muscle strength is more important than the skills for the weight lifting sports. And it is interesting to notice that the property of the variable stiffness actuators has to meet the core requirement in the different applications, like stiffness variation for the explosive movement or periodic movement and larger torque saturation in the same size and input power for the weight lifting movement.


## 6. Conclusion

In this article, the weight lifting problem of variable stiffness actuated robot was formulated as trajectory optimization problem. The original problem is discretized with Gauss pseudospectral method and solved with general sparse solver. Compared with the iLQR approach, the obtained trajectory is more smooth and highly precise. In addition, optimizations results with different setup of weight factors were presented to lead the choice of different optimal control policy. At last, the influence of stiffness variation on the weight lifting motion was discussed.In the future work, we intend to investigate the optimal control method through the experiment on variable stiffness actuated robot platform.



---

*Source: 4872189-2019-07-18.xml*