Discretized Mid-Value CLVI-PDNN Based Redundancy Resolution for Single Leg of Quadruped Robot

Abstract

The two most important performance indicators of quadruped robot are load capacity and walking speed, and these performance indicators of the whole robot finally reflect on the joint torques and angular velocities. To satisfy different requirements of walking speed and load capacity when quadruped robots implement different tasks, the joint torques and angular velocities need to be balanced with physical constraints of the joints. A single leg with redundant DOF (degree of freedom) could optimize the distribution of joint torques or angular velocities based on different performance requirements. This paper presents a kind of new recurrent neural networks taking joint torques and angular velocities simultaneously into consideration and proposes mid-value CLVI-PDNN to achieve the optimal joint torques and angular velocities with physical constraints of the mechanism as described in our previous paper. Because the continuous mid-value CLVI-PDNN has difficulty in real-time operation because of too much calculation workload, two kinds of methods are proposed to discretize the mid-value CLVI-PDNN for application on computer or digital circuit. The simulation results demonstrate the efficacy of the algorithm proposed in this paper.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Legged robots could be applied for rescue or service in unstructured environment due to its agility on rough terrain. Because of much more motion stability than biped robots and less complexity compared to hexapod robots or octopod robots, quadruped robots have a broad application prospect. Figure 1 shows the hydraulic actuated quadruped robot developed by our laboratory, which has four legs assembled symmetrically, four revolute joints on each leg, a torso with computer, and other electronic equipment on board.

[figure omitted; refer to PDF]

The schematic diagram of the single leg is shown in Figure 2, which has a redundant joint in the sagittal plane. The moment arms of actuating forces are marked in Figure 2 using red lines. If the performance of the hydraulic actuator is finitude, the maximum joint torques and angular velocity are determined by the length of these red lines. Dimension design of the moment arms of actuating forces directly affects the performance of walking speed and load capacity of the robot.

[figure omitted; refer to PDF]

Designing the mechanism dimensions of the legs especially the moment arms needs to optimize redundancy resolution of the single leg depending on performance demands then distribute joint torques and angular velocity on the condition of joint physical constraints. Simultaneously, optimizing the distribution of joint torques and angular velocity could be used for designing the control algorithm of the robot. Generally, two kinds of methods are widely used to deal with inverse kinematics of redundant robots, i.e., analytical method and numerical analysis.

In some analytical methods, inverse kinematics of redundant robot could be described as a parameterized function with parameters including joint angles [1, 2], the angle between the reference plane and robot [3, 4]. Manually adding constraints is another analytical method to calculate the inverse kinematics of redundant robot [5]. The inverse kinematics problem could be analytically solved while the redundant DOFs matches number of added constraints.

Because redundant robot systems are generally strongly nonlinear and highly coupled, it is not a simple task to resolve inverse kinematics problem. Redundancy resolution is not only beneficial to obstacles avoidance and mechanism physical constraints avoidance [6–8] but also is used for improving motion performance based on diverse criteria [9–12]. Pseudoinverse type solution is a kind of widely used redundancy resolution [13, 14]. Some researchers solved inverse kinematics problem using different methods including interpolation [15, 16] and gradient projection [17–19], but these methods were restricted to computation speed and error. A number of novel or so called intelligent numerical resolving methods are used for redundant inverse kinematics problem, including quadratic programming (QP), artificial neural networks, quaternion method, online learning algorithm, and genetic algorithm [20–24]. Among these methods, quadratic programming is tractable and has well expansibility to real-time application. Minimizing sum of squares of joint velocity is the most widely utilized optimization criterion which is computationally tractable but ignores dynamics at all [25]. Besides, the optimization criterion that minimizes sum of squares of joint torques is also exploited to resolve redundant inverse kinematics which may lead to instability.

To acquire a comprehensive performance combining kinematics and dynamics, Ma proposed an optimization criterion simultaneously minimizing joint velocity and torques [26]. Yunong Zhang also proposed an optimization criterion that consists of minimizing sum of squares of joint velocity and torques simultaneously [27]. The optimization problem on different level (i.e., joint velocity and joint torques) was converted into a standard QP problem with equality and inequality constraints on joint-acceleration level. Differing from pseudoinverse type solution used by Ma, Zhang numerically solved the constrained QP problem using a kind of recurrent neural network which is named linear variational inequality-based primal-dual neural network (LVI-PDNN) [28–30].

It is different from manipulators that legged robots alter between support phase and flight (swing) phase when they are walking or running. Therefore, the legs need to swing forward and support body in different phases alternatively. Accordingly, the legs need to provide sufficient joint velocity and joint torques in different phases alternatively. On the basis of the above analysis, a new kind of optimization criterion consisting of support phase and flight (swing) phase is proposed in our previous paper and converted onto angular acceleration level based on ZNN. A feasible distribution of the joint torques and angular velocity is obtained according to continuous mid-value CLVI-PDNN, but the algorithm used in our previous paper has a drawback, i.e., hard to real-time application because of too large computation cost [31].

For applying recurrent neural networks based real-time optimization on computer or digital circuit, discrete neural networks gradually attract attention from researchers. M.J. Perez-Ilzarbe proposed a kind of discrete recurrent neural network to solve bounded strictly convex quadratic programming problem [32]. Then he applied Wolf duality theory to form a kind of simple discrete recurrent neural network for solving linear constraints quadratic programming problem [33]. H. Tang et al. proposed a general recurrent neural network to solve hybrid constraints quadratic programming problem [34]. Q. Liu et al. used discrete recurrent neural network to solve equality constraints quadratic programming problem and proved the global exponential stability of its solution [35]. Y. Zhang, B. Liao, L. Jin et al. researched on a series of problems such as time-variant matrix inversion, time-variant polynomial root searching, and time-variant nonlinear optimization. They compared several discretization methods and proved the convergency of discrete recurrent neural network [36–43].

The main contributions of this paper are as follows: firstly, two kinds of discretization methods are proposed to discretize the mid-value CLVI-PDNN for computer control of the quadruped robot, i.e., bilinear transform-type and Taylor-type discretization methods, while the recursion formulas of the two kinds of discrete mid-value CLVI-PDNN are also presented; secondly, inherent nature of the two kinds of discretization algorithms is analyzed theoretically, including the order of truncation error and stability; thirdly, simulation results of the proposed methods indicate the efficacy for redundancy resolution of single leg.

2. The Optimization Criterion and QP Formulation

In this paper, the redundant single leg of the quadruped robot is illustrated to show how to deal with its inverse kinematics problem using optimization method. Firstly, the kinematics and dynamics of the single leg system in sagittal plane with three joints are presented in this section. The kinematics of single leg is as follows:$$\begin{array}{}\text{(1)}& r=f\left(\theta \right),\end{array}$$where $\theta ={\left[\begin{array}{ccc}{\theta }_{\mathrm{1}}& {\theta }_{\mathrm{2}}& {\theta }_{\mathrm{3}}\end{array}\right]}^{\mathrm{T}}$ is the position of joints in joint space, while $r={\left[\begin{array}{cc}{r}_{\mathrm{1}}& r_{\mathrm{2}}\end{array}\right]}^{\mathrm{T}}$ is the position of foot in Cartesian space. The differential kinematics is as follows:$$\begin{array}{}\text{(2)}& \dot{r}=J\left(\theta \right)\dot{\theta },\end{array}$$where $\dot{\theta }={\left[\begin{array}{ccc}{\dot{\theta }}_{\mathrm{1}}& {\dot{\theta }}_{\mathrm{2}}& {\dot{\theta }}_{\mathrm{3}}\end{array}\right]}^{\mathrm{T}}$ is the velocity of joints in joint space, while $\dot{r}={\left[\begin{array}{cc}{\dot{r}}_{\mathrm{1}}& {\dot{r}}_{\mathrm{2}}\end{array}\right]}^{\mathrm{T}}$ is the velocity of foot in Cartesian space. $J=\partial f/\partial \theta$ is Jacobian matrix. Second order differential equation is $$\begin{array}{}\text{(3)}& \ddot{r}=J\left(\theta \right)\ddot{\theta }+\dot{J}\left(\theta \right)\dot{\theta },\end{array}$$where $\ddot{\theta }={\left[\begin{array}{ccc}{\ddot{\theta }}_{\mathrm{1}}& {\ddot{\theta }}_{\mathrm{2}}& {\ddot{\theta }}_{\mathrm{3}}\end{array}\right]}^{\mathrm{T}}$ is the acceleration of joint in joint space, while $\ddot{r}={\left[\begin{array}{cc}{\ddot{r}}_{\mathrm{1}}& {\ddot{r}}_{\mathrm{2}}\end{array}\right]}^{\mathrm{T}}$ is the acceleration of foot in Cartesian space. $\dot{J}=\partial J/\partial t$ is differentiation of J with respect to time t. To simplify the description, Equation (3) is rewritten as follows:$$\begin{array}{}\text{(4)}& J\left(\theta \right)\ddot{\theta }={\ddot{r}}_a,\end{array}$$in which ${\ddot{r}}_a=\ddot{r}-\dot{J}\dot{\theta }$.

As the general form of robot dynamics, the single leg also is presented as follows:$$\begin{array}{}\text{(5)}& M\left(\theta \right)\ddot{\theta }+C\left(\theta ,\dot{\theta }\right)+G\left(\theta \right)=\tau ,\end{array}$$in which $M(\theta )$, $C(\theta ,\dot{\theta })$, and $G(\theta )$ denote inertia matrix, centripetal force and Coriolis force, and the gravitational term, respectively. $\tau ={\left[\begin{array}{ccc}{\tau }_{\mathrm{1}}& {\tau }_{\mathrm{2}}& {\tau }_{\mathrm{3}}\end{array}\right]}^{\mathrm{T}}$ is joint torques of the single leg.

It has been analyzed that a single optimization criterion is not suitable for legged robots that alter between support phase and flight (or swing) phase when walking or running. Then a modified optimization criterion that minimizes two norms including joint velocity and torques was proposed in our previous paper as follows:$$\begin{array}{}\text{(6)}& \mathrm{m}\mathrm{i}\mathrm{n}\left\{\begin{array}{cc}\frac{\mathrm{1}}{\mathrm{2}}\left({\alpha }_s{\dot{\theta }}^{\mathrm{T}}V\dot{\theta }+{\beta }_s{\tau }^{\mathrm{T}}W\tau \right),\hfill & \mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{t}\hspace{0.17em}\hspace{0.17em}\mathrm{p}\mathrm{h}\mathrm{a}\mathrm{s}\mathrm{e}\hfill \\ \frac{\mathrm{1}}{\mathrm{2}}\left({\alpha }_f{\dot{\theta }}^{\mathrm{T}}V\dot{\theta }+{\beta }_f{\tau }^{\mathrm{T}}W\tau \right),\hfill & \mathrm{f}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}\hspace{0.17em}\hspace{0.17em}\mathrm{p}\mathrm{h}\mathrm{a}\mathrm{s}\mathrm{e},\hfill \end{array}\right.\end{array}$$with equality constraints Equations (4) and (5) and inequality constraints, i.e., physical constraints of the joint as follows:$$\begin{array}{c}\text{(7)}& {\theta }^{-}\le \theta \le {\theta }^{+}{\dot{\theta }}^{-}\le \dot{\theta }\le {\dot{\theta }}^{+}\\ {\ddot{\theta }}^{-}\le \ddot{\theta }\le {\ddot{\theta }}^{+},\end{array}$$${\alpha }_s$ and ${\beta }_s$ are weight factors in support phase, while ${\alpha }_f$ and ${\beta }_f$ are weight factors in flight (or swing) phase. The values of these weight factors fall in the interval $\left[\begin{array}{cc}\mathrm{0}& \mathrm{1}\end{array}\right]$ and satisfy ${\alpha }_s+{\beta }_s=\mathrm{1}$, ${\alpha }_f+{\beta }_f=\mathrm{1}$. $V$ and $W$ are diagonal weight matrices reflecting effects of different joints.

The optimization problem formed by Equations (4)-(7) is converted into a standard QP problem on the joint angular acceleration level using the method base on ZNN in [27] as follows:$$\begin{array}{}\text{(8)}& \mathrm{m}\mathrm{i}\mathrm{n}\left\{\begin{array}{cc}\frac{\mathrm{1}}{\mathrm{2}}{\ddot{\theta }}^{\mathrm{T}}{Q}_s\ddot{\theta }+p_s^{\mathrm{T}}\ddot{\theta },\hfill & \mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{t}\hspace{0.17em}\hspace{0.17em}\mathrm{p}\mathrm{h}\mathrm{a}\mathrm{s}\mathrm{e}\hfill \\ \frac{\mathrm{1}}{\mathrm{2}}{\ddot{\theta }}^{\mathrm{T}}{Q}_f\ddot{\theta }+p_f^{\mathrm{T}}\ddot{\theta },\hfill & \mathrm{f}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}\hspace{0.17em}\hspace{0.17em}\mathrm{p}\mathrm{h}\mathrm{a}\mathrm{s}\mathrm{e},\hfill \end{array}\right.\end{array}$$subjected to equality and inequality constraints$$\begin{array}{}\text{(9)}& {\ddot{r}}_a=J\ddot{\theta },\hspace{1em}{\ddot{\theta }}^{-}\le \ddot{\theta }\le {\ddot{\theta }}^{+},\end{array}$$in which$$\begin{array}{c}\text{(10)}& Q_s={\alpha }_{s}V+\left(\mathrm{1}-{\alpha }_s\right)M^{\mathrm{T}}WM{p}_s^{\mathrm{T}}={\alpha }_s\lambda V\dot{\theta }+\left(\mathrm{1}-{\alpha }_s\right)M^{\mathrm{T}}{W}^{\mathrm{T}}\left(C+G\right),\\ \text{(11)}& Q_f={\alpha }_{f}V+\left(\mathrm{1}-{\alpha }_f\right)M^{\mathrm{T}}WM{p}_f^{\mathrm{T}}={\alpha }_f\lambda V\dot{\theta }+\left(\mathrm{1}-{\alpha }_f\right)M^{\mathrm{T}}{W}^{\mathrm{T}}\left(C+G\right),\end{array}$$$\lambda$ is a positive constant corresponding to convergence rate of ZNN.

The standard QP problem formulated by Equations (8) and (9) is equivalent to a kind of neural network named LVI-PDNN (LVI-based primal-dual neural network) as follows:$$\begin{array}{}\text{(12)}& \dot{y}=\left\{\begin{array}{cc}\eta \left(I+H_s^{\mathrm{T}}\right)\left\{P_{\Omega }\left[y-\left(H_{s}y+q_s\right)\right]-y\right\},\hfill & \mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{t}\hspace{0.17em}\hspace{0.17em}\mathrm{p}\mathrm{h}\mathrm{a}\mathrm{s}\mathrm{e}\hfill \\ \eta \left(I+H_f^{\mathrm{T}}\right)\left\{P_{\Omega }\left[y-\left(H_{f}y+q_f\right)\right]-y\right\},\hfill & \mathrm{f}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}\hspace{0.17em}\hspace{0.17em}\mathrm{p}\mathrm{h}\mathrm{a}\mathrm{s}\mathrm{e},\hfill \end{array}\right.\end{array}$$where $y=[{\ddot{\theta }}^{\mathrm{T}},u^{\mathrm{T}}{]}^{\mathrm{T}},\hspace{0.17em}{q}_s=[p_s,-{\ddot{r}}_a^{\mathrm{T}}{]}^{\mathrm{T}},\hspace{0.17em}{q}_f=[p_f,-{\ddot{r}}_a^{\mathrm{T}}{]}^{\mathrm{T}},\hspace{0.17em}{H}_s=\left[\begin{array}{cc}{Q}_s& -J_L^{\mathrm{T}}\\ J_L& \mathrm{0}\end{array}\right],\hspace{0.17em}{H}_f=\left[\begin{array}{cc}{Q}_f& -J_L^{\mathrm{T}}\\ J_L& \mathrm{0}\end{array}\right]$, and $u$ denotes the dual decision variable vector consisting of coefficients of equality constraints in Equation (9). $\eta$ is a positive constant that determines convergence rate of LVI-PDNN and should be selected large enough. $P_{\Omega }(·)$ is a piecewise linear projection operator as follows:$$\begin{array}{c}\text{(13)}& P_{\Omega }\left(·\right):R^{n+m+l}\to \Omega =\left\{y\mid {\varsigma }^{-}\le y\le {\varsigma }^{+}\right\}{\left[P_{\Omega }\left(y\right)\right]}_i=\left\{\begin{array}{cc}{\varsigma }_i^{-},\hfill & y_i<{\varsigma }_i^{-}\hfill \\ y_i,\hfill & {\varsigma }_i^{-}\le y_i\le {\varsigma }_i^{+}\hfill \\ {\varsigma }_i^{+},\hfill & y_i>{\varsigma }_i^{+},\hfill \end{array}\right.\hspace{1em}i=\mathrm{1},\cdots ,n+m+l.\end{array}$$

3. CLVI-PDNN and Mid-Value CLVI-PDNN

In our previous paper [31] two kinds of neural network are proposed based on LVI-PDNN, i.e., CLVI-PDNN and mid-value CLVI-PDNN. CLVI-PDNN improves the drawback of open-loop LVI-PDNN such as divergency of the optimized results from the desired value. Mid-value CLVI-PDNN further solves the problem that the optimized joint angles violate the physical constraints of the mechanism and deviate from the desired configuration especially for a legged robot.

CLVI-PDNN means closed-loop LVI-PDNN that introduces the output to the input of LVI-PDNN as feedback shown in Figure 3. In Figure 3, the block diagram in the blue dotted line shows LVI-PDNN proposed by Zhang [27]. The block diagram outside of the blue dotted line is a closed loop introducing feedback to the input of the neural network system as follows:$$\begin{array}{}\text{(14)}& {\ddot{r}}_a={\ddot{r}}_d+k_{pr}\left(r_d-f\left(\theta \right)\right)+k_{ir}{{\displaystyle \int }}_{\mathrm{0}}^t\left(r_d-f\left(\theta \right)\right)dt+k_{pv}\left({\dot{r}}_d-J\left(\theta \right)\dot{\theta }\right)+k_{iv}{{\displaystyle \int }}_{\mathrm{0}}^t\left({\dot{r}}_d-J\left(\theta \right)\dot{\theta }\right)dt-\dot{J}\left(\theta \right)\dot{\theta },\end{array}$$in which $r_d$, ${\dot{r}}_d$, and ${\ddot{r}}_d$ denote the desired motion of the foot in Cartesian space, respectively, while $k_{pr}$, $k_{ir}$, $k_{pv}$, and $k_{iv}$ are coefficient matrices.

[figure omitted; refer to PDF]

CLVI-PDNN has a critical defect that the optimized joint angles exceed the joint physical limitations leading to its uselessness for motion control of the quadruped robot. To improve this weakness of CLVI-PDNN, it could be modified as follows: defining a slack variable $\epsilon =\theta -{\theta }_d$, where ${\theta }_d$ denotes the desired joint configuration of the single leg when walking or running. $\epsilon$ could keep stable using a simple feedback control law$$\begin{array}{}\text{(15)}& \theta ={\theta }_{\mathrm{0}}+k_p\left({\theta }_d-\theta \right)+{{\displaystyle \int }}_{\mathrm{0}}^t\dot{\theta }dt,\end{array}$$in which ${\theta }_{\mathrm{0}}$ means initial values and $k_p$ means positive coefficient matrix. To prevent joint angles from exceeding the joint physical constraints, value of ${\theta }_d$ should be selected close to the mid-value of the joint motion range. That is also why this algorithm is named mid-value CLVI-PDNN. Figure 4 shows the block diagram of mid-value CLVI-PDNN.

[figure omitted; refer to PDF]

As analyzed in our previous paper [31], the mid-value CLVI-PDNN just compromise between the inequality and equality constrains. That is to say, the output of the QP solver maintains in motion range of the joints ensuring a feasible rather than optimal solution on condition that the error between the desired and the optimized trajectory is not enough to destroy normal motion control. To some extent, in many situations especially when the quadruped robots walk or run in the field, the output of mid-value CLVI-PDNN would be useful for motion control.

For some circumstances, such as high dimensions of the augmented matrix $M_s$ and $M_f$, restricted performance of the onboard controller or computer, the discrete form of the QP solvers mid-value CLVI-PDNN may be more applicable for real-time control. We will deduce the discrete form of mid-value CLVI-PDNN and analyze the characteristics in next section.

4. Discretize Mid-Value CLVI-PDNN

To solve the problem that continuous mid-value CLVI-PDNN is difficult to be applied in real-time control system, two kinds of discretization methods are presented in this section. For a general unconstraint optimization problem$$\begin{array}{}\text{(16)}& \underset{x\left(t\right)\in R^n}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}f\left(x\left(t\right),t\right).\end{array}$$Suppose the extreme value of equation as follows:$$\begin{array}{}\text{(17)}& {\Omega }_{\ast }\left(t\right)=\left\{\left(x^{\ast }\left(t\right),t\right)\mid g\left(x^{\ast }\left(t\right),t\right)=\frac{\partial f\left(x^{\ast }\left(t\right),t\right)}{\partial x^{\ast }\left(t\right)}=\mathrm{0}\right\}.\end{array}$$Obviously, condition of the extreme value is $g(x(t),t)=\mathrm{0}$, so that the error equation is selected as $e(t)=g(x(t),t)$. Then the error dynamics is designed on the basis of ZNN as follows:$$\begin{array}{}\text{(18)}& \dot{e}\left(t\right)=-\lambda e\left(t\right),\end{array}$$where $\lambda$ is a positive constant to scale the rate of convergency of $e(t)$. Substituting $e(t)=g(x(t),t)$ into Equation (18) could obtain$$\begin{array}{}\text{(19)}& H\left(x\left(t\right),t\right)\dot{x}\left(t\right)=-\lambda g\left(x\left(t\right),t\right)-g^{\prime }\left(x\left(t\right),t\right),\end{array}$$where$$\begin{array}{c}\text{(20)}& H\left(x\left(t\right),t\right)=\frac{\partial g\left(x\left(t\right),t\right)}{\partial x^T}=\frac{{\partial }^{\mathrm{2}}f\left(x\left(t\right),t\right)}{\partial x\left(t\right)\partial x^T\left(t\right)},g^{\prime }\left(x\left(t\right),t\right)=\frac{\partial g\left(x\left(t\right),t\right)}{\partial t}=\frac{{\partial }^{\mathrm{2}}f\left(x\left(t\right),t\right)}{\partial x\left(t\right)\partial t}.\end{array}$$If the Hessian matrix $H(x(t),t)$ is nonsingular, multiplying $H^{-\mathrm{1}}(x(t),t)$ both sides of Equation (19) results in$$\begin{array}{}\text{(21)}& \dot{x}\left(t\right)=-H^{-\mathrm{1}}\left(x\left(t\right),t\right)\left(\lambda g\left(x\left(t\right),t\right)+g^{\prime }\left(x\left(t\right),t\right)\right).\end{array}$$Given initial value $x(\mathrm{0})$, the optimal state $x^{\ast }(t)$ will be achieved on condition that Hessian matrix $H(x(t),t)$ is positive definite. Actually, for our quadruped robot, Hessian matrix always maintains positive definite.

Taking forward Euler method as example, the discretization process could be presented as follows: rewriting Equation (21) as $\dot{x}=G(x(t),t)$, the forward Euler equation is as follows:$$\begin{array}{}\text{(22)}& x_{k+\mathrm{1}}=x_k+\delta G\left(x_k,t_k\right),\end{array}$$where $x_k=x(t_k)$ means the value of $x(t)$ at moment $t_k$. $t_{k+\mathrm{1}}=t_k+\delta$, in which $\delta$ denotes the sample interval. Therefore, discretization equation of $x(t)$ is as follows:$$\begin{array}{}\text{(23)}& x_{k+\mathrm{1}}=x_k-H^{-\mathrm{1}}\left(x_k,t_k\right)\left(\delta \lambda g\left(x_k,t_k\right)+\delta g^{\prime }\left(x_k,t_k\right)\right).\end{array}$$Next the above-mentioned discretization process is applied for redundancy resolution of single leg. For a time-variant QP problem, it could be described at moment $t_k$ as follows:$$\begin{array}{c}\text{(24)}& \mathrm{m}\mathrm{i}\mathrm{n}\hspace{1em}\frac{\mathrm{1}}{\mathrm{2}}{\ddot{\theta }}_k^T{Q}_k{\ddot{\theta }}_k+p_k^T{\ddot{\theta }}_k\mathrm{s}.\mathrm{t}.\hspace{1em}{J}_L^k{\ddot{\theta }}_k={\ddot{r}}_a^k\\ \hspace{1em}{\xi }^{-}\le {\ddot{\theta }}_k\le {\xi }^{+}.\end{array}$$Moreover, the dynamics of LVI-PDNN at moment $t_k$ is as follows:$$\begin{array}{}\text{(25)}& {\dot{y}}_k=\eta \left(I+H_k^{\mathrm{T}}\right)\left\{P_{\Omega }\left[y_k-\left(H_k{y}_k+q_k\right)\right]-y_k\right\},\end{array}$$where $y_k=\left[\begin{array}{c}{\ddot{\theta }}_k\\ u_k\end{array}\right],\hspace{0.17em}{q}_k=\left[\begin{array}{c}{p}_k\\ -{\ddot{r}}_a^k\end{array}\right],\hspace{0.17em}{H}_k=\left[\begin{array}{cc}{Q}_k& -J_L^{k\mathrm{T}}\\ J_L^k& \mathrm{0}\end{array}\right]$.

For different discretization methods, accuracy and rapidity differ from each other. Taylor-type and bilinear transform-type discretization methods are used for mid-value CLVI-PDNN in this paper. The discretized state of neural network using bilinear transform is as follows:$$\begin{array}{}\text{(26)}& y_{k+\mathrm{1}}=y_k+\frac{\delta }{\mathrm{2}}\left({\dot{y}}_{k+\mathrm{1}}+{\dot{y}}_k\right).\end{array}$$Substituting Equation (25) into Equation (26) obtains$$\begin{array}{}\text{(27)}& y_{k+\mathrm{1}}=y_k+\frac{\delta }{\mathrm{2}}\eta \left(I+H^{\mathrm{T}}\right)\left(\left\{P_{\Omega }\left[y_{k+\mathrm{1}}-\left(H_{k+\mathrm{1}}{y}_{k+\mathrm{1}}+q_{k+\mathrm{1}}\right)\right]-y_{k+\mathrm{1}}\right\}+\left\{P_{\Omega }\left[y_k-\left(H_k{y}_k+q_k\right)\right]-y_k\right\}\right).\end{array}$$Then, the recursion formula of $y_k$ could be achieved:$$\begin{array}{}\text{(28)}& y_{k+\mathrm{1}}=\frac{\mathrm{2}-\delta }{\mathrm{2}+\delta }{y}_k+\frac{\delta }{\mathrm{2}+\delta }\eta \left(I+H^{\mathrm{T}}\right)\left(P_{\Omega }\left[y_{k+\mathrm{1}}-\left(H_{k+\mathrm{1}}{y}_{k+\mathrm{1}}+q_{k+\mathrm{1}}\right)\right]+P_{\Omega }\left[y_k-\left(H_k{y}_k+q_k\right)\right]\right).\end{array}$$$y_{k+\mathrm{1}}$ in the right-hand side of Equation (28) could be estimated according to forward difference$$\begin{array}{c}\text{(29)}& y_{k+\mathrm{1}}=y_k+\delta {\dot{y}}_k=y_k+\delta \eta \left(I+H^{\mathrm{T}}\right)\left\{P_{\Omega }\left[y_k-\left(H_k{y}_k+q_k\right)\right]-y_k\right\},H_{k+\mathrm{1}}=H_k+\delta {\dot{H}}_k=H_k+\Delta H_k,\\ q_{k+\mathrm{1}}=q_k+\delta {\dot{q}}_k=q_k+\Delta q_k,\end{array}$$where $\Delta H_k=H_k-H_{k-\mathrm{1}}$ and $\Delta q=q_k-q_{k-\mathrm{1}}$.

Taylor-type discretization method exploits high order terms of the Taylor series and then utilizes multiple points to discretize mid-value CLVI-PDNN. This kind of discretization method possesses the advantage of higher computational accuracy. Higher order corresponds higher accuracy while consuming more computation. Eliminating third-order derivative terms of Taylor series obtains$$\begin{array}{}\text{(30)}& {\dot{y}}_k=\frac{\mathrm{2}{y}_{k+\mathrm{1}}-\mathrm{3}{y}_k+\mathrm{2}{y}_{k-\mathrm{1}}-y_{k-\mathrm{2}}}{\mathrm{2}\delta }+O\left({\delta }^{\mathrm{2}}\right).\end{array}$$Substituting Equation (25) into Equation (30) achieves the recursion formula as follows:$$\begin{array}{}\text{(31)}& y_{k+\mathrm{1}}=\mathrm{1.5}{y}_k-y_{k-\mathrm{1}}+\mathrm{0.5}{y}_{k-\mathrm{2}}+\delta {\dot{y}}_k=\mathrm{1.5}{y}_k-y_{k-\mathrm{1}}+\mathrm{0.5}{y}_{k-\mathrm{2}}+\delta \eta \left(I+H_k^{\mathrm{T}}\right)\left\{P_{\Omega }\left[y_k-\left(H_k{y}_k+q_k\right)\right]-y_k\right\}.\end{array}$$

5. Theory Analysis

In this section, convergency and steady-state residual error of the discretization methods used in this paper is analyzed.

5.1. Bilinear Transform-Type Discretization Method

Theorem 1.

Bilinear transform-type recursion formula as Equation (28) is 0-stable.

Proof.

The characteristic polynomial of bilinear transform-type recursion formula as Equation (28) is$$\begin{array}{}\text{(32)}& P\left(z\right)=z-\frac{\mathrm{2}-\delta }{\mathrm{2}+\delta }=\mathrm{0},\end{array}$$which has a unique root on/in the unit disk, i.e., $z=(\mathrm{2}-\delta )/(\mathrm{2}+\delta )$ with $\mathrm{0}<\delta <\mathrm{1}$. Therefore, bilinear transform-type recursion formula as Equation (28) is 0-stable. The proof is thus completed.

Theorem 2.

Bilinear transform-type recursion formula as Equation (28) is consistent and convergent, which converges with the order of truncation error $O({\delta }^{\mathrm{2}})$.

Proof.

Discretization formula of $y_k$ based on forward difference formula is as follows:$$\begin{array}{}\text{(33)}& {\dot{y}}_k=\frac{y_{k+\mathrm{1}}-y_k}{\delta }+O\left(\delta \right).\end{array}$$While discretization formula of $y_k$ based on backward difference formula is as follows:$$\begin{array}{}\text{(34)}& {\dot{y}}_{k+\mathrm{1}}=\frac{y_{k+\mathrm{1}}-y_k}{\delta }+O\left(\delta \right).\end{array}$$Adding Equations (33) and (34) could obtain the bilinear transform-type discretization formula$$\begin{array}{}\text{(35)}& y_{k+\mathrm{1}}=y_k+\frac{\delta }{\mathrm{2}}\left({\dot{y}}_{k+\mathrm{1}}+{\dot{y}}_k\right)+O\left({\delta }^{\mathrm{2}}\right).\end{array}$$Note that dropping $O({\delta }^{\mathrm{2}})$ from Equation (35) yields exactly bilinear transform-type recursion formula as Equation (26), and, thus, the truncation error of Equation (26) is $O({\delta }^{\mathrm{2}})$. Consequently, bilinear transform-type recursion formula as Equation (28) is consistent with order 2. Combining with Theorem 1, we know that the model in Equation (28) is both 0-stable and consistent, so that it can be conclusively concluded that bilinear transform-type recursion formula as Equation (28) is consistent and convergent, which converges with the order of its truncation error being $O({\delta }^{\mathrm{2}})$. The proof is now completed.

Theorem 3.

Considering the QP solver mid-value CLVI-PDNN based on ZNN, the steady-state residual error ${\lim }_{k\to \infty }{‖{\dot{y}}_k‖}_{\mathrm{2}}$ of bilinear transform-type discretization recursion formula as Equation (28) is $O({\delta }^{\mathrm{2}})$ in the sense of Euclidean norm.

Proof.

Let $y_k^{\ast }$ denote the theoretical optimal solution of the QP problem described by Equations (8) and (9), we know that $P_{\Omega }[y_k^{\ast }-(H_k{y}_k^{\ast }+q_k)]-y_k^{\ast }=\mathrm{0}$. On the basis of Theorems 1 and 2, it could be concluded that $y_k=y_k^{\ast }+O({\delta }^{\mathrm{2}})$ with $k$ large enough. Therefore, we obtain the limit error (the definition of error follows [42]):$$\begin{array}{}\text{(36)}& \underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖{\dot{y}}_k‖}_{\mathrm{2}}=\underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖\eta \left(I+H_k^{\mathrm{T}}\right)\left\{P_{\Omega }\left[y_k-\left(H_k{y}_k+q_k\right)\right]-y_k\right\}‖}_{\mathrm{2}}=\underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖\eta \left(I+H_k^{\mathrm{T}}\right)\left\{P_{\Omega }\left[y_k^{\ast }+O\left({\delta }^{\mathrm{2}}\right)-\left(H_k\left(y_k^{\ast }+O\left({\delta }^{\mathrm{2}}\right)\right)+q_k\right)\right]-y_k^{\ast }-O\left({\delta }^{\mathrm{2}}\right)\right\}‖}_{\mathrm{2}}=\underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖\eta \left(I+H_k^{\mathrm{T}}\right)\left\{P_{\Omega }\left[y_k^{\ast }-\left(H_k{y}_k^{\ast }+q_k\right)\right]-y_k^{\ast }\right\}+\left(\eta \left(I+H_k^{\mathrm{T}}\right)\left(I-H_k\right)-I\right)O\left({\delta }^{\mathrm{2}}\right)‖}_{\mathrm{2}}=\underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖\left(\eta \left(I+H_k^{\mathrm{T}}\right)\left(I-H_k\right)-I\right)O\left({\delta }^{\mathrm{2}}\right)‖}_{\mathrm{2}}.\end{array}$$Considering that $((I+H_k^{\mathrm{T}})(I-H_k)-I)$ is a constant matrix at a certain time instant, we further obtain$$\begin{array}{}\text{(37)}& \underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖{\dot{y}}_k‖}_{\mathrm{2}}=\underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖\left(\eta \left(I+H_k^{\mathrm{T}}\right)\left(I-H_k\right)-I\right)O\left({\delta }^{\mathrm{2}}\right)‖}_{\mathrm{2}}=O\left({\delta }^{\mathrm{2}}\right).\end{array}$$The proof is completed.

5.2. Taylor-Type Discretization Method

Taylor-type recursion formula as Equation (30) is 0-stable. Taylor-type recursion formula as Equation (30) is consistent and convergent, which converges with the order of truncation error $O({\delta }^{\mathrm{3}})$. These conclusions are proved in [42]. Therefore, the steady-state residual error could be concluded as follows.

Theorem 4.

Considering the QP solver mid-value CLVI-PDNN based on ZNN, the steady-state residual error ${\lim }_{k\to \infty }{‖{\dot{y}}_k‖}_{\mathrm{2}}$ of Taylor-type discretization recursion formula as Equation (31) is $O({\delta }^{\mathrm{3}})$ in the sense of Euclidean norm.

Proof.

Let $y_k^{\ast }$ denote the theoretical optimal solution of the QP problem described by Equations (8) and (9), so that $P_{\Omega }[y_k^{\ast }-(H_k{y}_k^{\ast }+q_k)]-y_k^{\ast }=\mathrm{0}$ and $y_k=y_k^{\ast }+O({\delta }^{\mathrm{3}})$ with k large enough. Therefore, we obtain the limit error as follows:$$\begin{array}{}\text{(38)}& \underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖{\dot{y}}_k‖}_{\mathrm{2}}=\underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖\eta \left(I+H_k^{\mathrm{T}}\right)\left\{P_{\Omega }\left[y_k-\left(H_k{y}_k+q_k\right)\right]-y_k\right\}‖}_{\mathrm{2}}=\underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖\eta \left(I+H_k^{\mathrm{T}}\right)\left\{P_{\Omega }\left[y_k^{\ast }+O\left({\delta }^{\mathrm{3}}\right)-\left(H_k\left(y_k^{\ast }+O\left({\delta }^{\mathrm{3}}\right)\right)+q_k\right)\right]-y_k^{\ast }-O\left({\delta }^{\mathrm{3}}\right)\right\}‖}_{\mathrm{2}}=\underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖\eta \left(I+H_k^{\mathrm{T}}\right)\left\{P_{\Omega }\left[y_k^{\ast }-\left(H_k{y}_k^{\ast }+q_k\right)\right]-y_k^{\ast }\right\}+\left(\eta \left(I+H_k^{\mathrm{T}}\right)\left(I-H_k\right)-I\right)O\left({\delta }^{\mathrm{3}}\right)‖}_{\mathrm{2}}=\underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}{‖\left(\eta \left(I+H_k^{\mathrm{T}}\right)\left(I-H_k\right)-I\right)O\left({\delta }^{\mathrm{3}}\right)‖}_{\mathrm{2}}=O\left({\delta }^{\mathrm{3}}\right).\end{array}$$The proof is completed.

6. Results and Discussion

The proposed discrete QP algorithms are operated in Matlab to solve the optimization problem formulated by Equations (8) and (9). In the simulation environment, the quadruped robot moves along horizontal direction with trotting gait in sagittal plane. The body of the quadruped robot moves at a constant speed on horizontal direction as well as constant height. The parameters of trotting gait are listed in Table 1. h is a constant meaning the height of the hip. $h_{apex}$ is the maximum height of the foot when the quadruped robot is walking. $v$ is a constant meaning the desired velocity of the body. $T$ is the gait cycle period of trotting gait. $g$ is the gravitational acceleration. Parameters of the single leg system with half mass of the whole body in the sagittal plane as shown in Figure 2 are listed in Table 2. $M$ is mass of body of the single system, i.e., half mass of the whole body. The others mean the masses and lengths of the three links of the single leg in sagittal plane.

Table 1

Parameters of trotting gait.

parameter	value	unit
$$h$$	0.7	m
$$h_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{x}}$$	0.15	m
$$v$$	6.0	Km/h
$$T$$	0.6	s

Table 2

Physical parameters of single leg.

parameter	value	unit
$$M$$	50	kg
$$M_1$$	3	kg
$$M_2$$	3	kg
$$M_3$$	3	kg
$$L_1$$	0.25	m
$$L_2$$	0.30	m
$$L_3$$	0.30	m

The joint physical constraints, i.e., the upper and lower bound of joint angles, angular velocity and acceleration are listed in Table 3. The parameters of the discretized mid-value CLVI-PDNN are listed in Table 4.

Table 3

Physical constraints of joints.

parameter	value	unit
$${\theta }^{+}$$	$\left(17\pi /36,5\pi /6,5\pi /6\right)$	rad
$${\theta }^{-}$$	$\left(\pi /9,\pi /9,\pi /9\right)$	rad
$${\dot{\theta }}^{+}$$	$\left(12,12,20\right)$	rad/s
$${\dot{\theta }}^{-}$$	$-\left(12,12,20\right)$	rad/s
$${\ddot{\theta }}^{+}$$	$$\left(1700,2200,2800\right)$$	rad/s
$${\ddot{\theta }}^{-}$$	$$-\left(1700,2200,2800\right)$$	rad/s

Table 4

Parameters of mid-value CLVI-PDNN.

parameter	value
$$\lambda$$	3600
$$\eta$$	107
$$k_{pr}$$	1000
$$k_{ir}$$	100000
$$k_{pv}$$	100000
$$k_{iv}$$	100000

In addition, motion trajectory of foot is designed as a curve respect to time $t$ on the horizontal and vertical direction respectively. The expected motion curves in horizontal direction and vertical direction are described as follows:$$\begin{array}{}\text{(39)}& x=\left\{\begin{array}{cc}-\frac{\pi }{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\mathrm{20}\pi }{\mathrm{3}}t\right)+\frac{\mathrm{5}}{\mathrm{3}}t-\frac{\mathrm{1}}{\mathrm{4}},\hfill & \mathrm{0}\le t\le \mathrm{0.3}\hfill \\ \frac{\mathrm{1}}{\mathrm{4}}-\frac{\mathrm{5}}{\mathrm{3}}\left(t-\mathrm{0.3}\right),\hfill & \mathrm{0.3}<t\le \mathrm{0.6},\hfill \end{array}\right.\text{(40)}& y=\left\{\begin{array}{cc}-\mathrm{0.7}+\mathrm{237}{t}^{\mathrm{3}}-\mathrm{2370.4}{t}^{\mathrm{4}}+\mathrm{7901.2}{t}^{\mathrm{5}}-\mathrm{8779.1}{t}^{\mathrm{6}},\hfill & \mathrm{0}\le t\le \mathrm{0.3}\hfill \\ -\mathrm{0.7},\hfill & \mathrm{0.3}<t\le \mathrm{0.6}.\hfill \end{array}\right.\end{array}$$The step length is set the same as the control cycle, i.e., $\delta =\mathrm{1}$ms. The simulation results of the single leg under the aforementioned conditions are illustrated in Figures 5–10, in which Figures 5–7 present simulation results acquired by using Taylor-type method to discretize the mid-value CLVI-PDNN, while Figures 8–10 present simulation results acquired with bilinear transform-type discretization method.

[figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF]

The parameter set 1 in Figure 5 are $V=W=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathrm{1,1},\mathrm{1}])$ ($\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(X)$ means a diagonal matrix whose diagonal elements are vector $X$), ${\alpha }_s=\mathrm{0.9},\hspace{0.17em}{\alpha }_f=\mathrm{0.1},\hspace{0.17em}{\theta }_d=(\mathrm{0.6379,1.8059,1.0472})$; the parameter set 2 in Figure 6 are $V=W=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathrm{1,1},\mathrm{1}]),\hspace{0.17em}{\alpha }_s=\mathrm{0.1},\hspace{0.17em}{\alpha }_f=\mathrm{0.9}$, ${\theta }_d=(\mathrm{0.6379,1.8059,1.0472})$; the parameter set 3 in Figure 7 are $V=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathrm{4,2},\mathrm{1}]),\hspace{0.17em}W=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathrm{1,2},\mathrm{4}])$, ${\alpha }_s=\mathrm{0.1}$, ${\alpha }_f=\mathrm{0.9}$, ${\theta }_d=\left(\begin{array}{ccc}\mathrm{0.6379}& \mathrm{1.8059}& \mathrm{1.0472}\end{array}\right)$. It is concluded that the optimized trajectory of the single leg in operation space could track the desired value well when using Taylor-type discretization method with position error about ±0.5cm in horizontal direction and ±2.5cm in vertical direction. For quadruped robot walking in the field that requires less position precision, the position error shown in Figures 5–7 could satisfy the design and control requirements relative to the motion range of the single leg. On the other hand, different parameter sets result in different optimized motion curves in joint space. As illustrated in Figures 5 and 6, the joint angular velocities in Figure 6 with ${\alpha }_s=\mathrm{0.1}$ are smoother than the curves in Figure 5 with ${\alpha }_s=\mathrm{0.9}$ demonstrating that norm of minimizing joint torques would lead to instability. The angular velocity of hip joint with $V=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathrm{1,1},\mathrm{1}])$ shown in Figure 6 is larger than the value with $V=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathrm{4,2},\mathrm{1}])$ shown in Figure 7. While angular velocity of ankle joint shown in Figure 6 is smaller than the value shown in Figure 7. Moreover, the motion ranges of joint angles in Figures 6 and 7 are quite distinct correspondingly. The comparison between Figures 6 and 7 indicates the effect of different velocity weights among the three joints.

The parameter set 1 in Figure 8 are $V=W=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathrm{1,1},\mathrm{1}])$, ${\alpha }_s=\mathrm{0.9},\hspace{0.17em}{\alpha }_f=\mathrm{0.1},\hspace{0.17em}{\theta }_d=(\mathrm{0.6379,1.8059,1.0472})$; the parameters set 2 in Figure 9 are $V=W=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathrm{1,1},\mathrm{1}]),\hspace{0.17em}{\alpha }_s=\mathrm{0.1}$, ${\alpha }_f=\mathrm{0.9}$, ${\theta }_d=(\mathrm{0.6379,1.8059,1.0472})$; the parameters set 3 in Figure 10 are $V=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathrm{4,2},\mathrm{1}]),\hspace{0.17em}W=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}([\mathrm{1,2},\mathrm{4}]),\hspace{0.17em}{\alpha }_s=\mathrm{0.1}$, ${\alpha }_f=\mathrm{0.9},\hspace{0.17em}{\theta }_d=\left(\begin{array}{ccc}\mathrm{0.6379}& \mathrm{1.8059}& \mathrm{1.0472}\end{array}\right)$. It is concluded that the optimized trajectory of the single leg in operation space could track the desired value well using bilinear transform-type discretization method with position error about ±0.5cm in horizontal direction and ±2.5cm in vertical direction. But it is disappointing that the motion trajectory in joint space is unsmooth and fluctuant so that perhaps it is difficult to be utilized in control of the single leg on the real system.

According to the simulation results of the two kinds of discretization methods, i.e., Taylor-type and bilinear transform-type, a conclusion could be achieved easily that the optimized motion trajectory in operation space tracks the desired value well and satisfies the control requirements. But the control of the single leg executes in joint space actually, the optimized motion trajectory in joint space by using Taylor-type discretization method is better than the other one obviously and easier to be applied in the control of the real mechanism.

7. Conclusion

As the discussion in our previous paper [31], we try to find a practicable QP solver to solve the inverse kinematics of the redundant single leg of the quadruped robot. In this paper we adopt two kinds of discretization methods to discretize the mid-value CLVI-PDNN proposed in our previous paper. According to the simulation results and analysis we could conclude that the optimized motion trajectories in operation space obtained by both Taylor-type and bilinear transform-type discretization methods track the desired values well if we do not overly concern about the position control precision of the single leg in the wild. Meanwhile, it is obvious that the optimized motion trajectory in joint space obtained by Taylor-type discretization method is superior to the other one. Of course, higher order Taylor-type discretization methods will achieve more accurate inverse kinematics of the redundant single leg.

It was analyzed in our previous paper that the mid-value CLVI-PDNN QP solver just outputs feasible solutions rather than intuitional optimal solutions corresponding to the optimization criterion formulated in Equation (6). In addition, the deviation between the optimized value and the desired value is much larger than theory analysis. The error may be caused by the constraints such as integration threshold, the deviation between the criteria in Equations (6) and (8), and probably the crucial reason, i.e., the introduction of expected joint configuration. Although the optimized value deviates from the desired value to some extent, it is satisfactory to achieve a useable answer in the end.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

This research was funded by the Preschool Research Program Funds for the National University of Defense Technology (No. ZK17-03-49).

References

[1] S. Lee, A. K. Bejczy, "Redundant arm kinematic control based on parameterization," Proceedings of the 1991 IEEE International Conference on Robotics and Automation, vol. 1, pp. 458-465, .

[2] D. Tolani, N. I. Badler, "Real-time inverse kinematics of the human arm," Presence: Teleoperators and Virtual Environments, vol. 5 no. 4, pp. 393-401, DOI: 10.1162/pres.1996.5.4.393, 1996.

[3] M. Shimizu, H. Kakuya, W.-K. Yoon, K. Kitagaki, K. Kosuge, "Analytical inverse kinematic computation for 7-DOF redundant manipulators with joint limits and its application to redundancy resolution," IEEE Transactions on Robotics, vol. 24 no. 5, pp. 1131-1142, DOI: 10.1109/TRO.2008.2003266, 2008.

[4] G. K. Singh, J. Claassens, "An analytical solution for the inverse kinematics of a redundant 7dof manipulator with link offsets," Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 2976-2982, .

[5] T. Zhang, Q. Wei, H. Ma, "Position/force control for a single leg of a quadruped robot in an operation space," International Journal of Advanced Robotic Systems, vol. 10 no. 137, 2013.

[6] H. Ding, S. P. Chan, "Real-time planning algorithm for obstacle avoidance of redundant robots," Journal of Intelligent & Robotic Systems, vol. 16 no. 3, pp. 229-243, 1996.

[7] A. Khoukhi, L. Baron, M. Balazinski, "A projected gradient augmented Lagrangian approach to multi-objective trajectory planning of redundant robots," Transactions of the Canadian Society for Mechanical Engineering, vol. 31 no. 4, pp. 391-405, DOI: 10.1139/tcsme-2007-0028, 2007.

[8] R. V. Mayorga, P. Sanongboon, "An artificial neural network approach for inverse kinematics computation and singularities prevention of redundant manipulators," Journal of Intelligent & Robotic Systems, vol. 44 no. 1, 2005.

[9] K. A. O'Neil, "Divergence of linear acceleration-based redundancy resolution schemes," IEEE Transactions on Robotics and Automation, vol. 18 no. 4, pp. 625-631, DOI: 10.1109/TRA.2002.801046, 2002.

[10] J. A. Goble, Y. Zhang, Y. Shimansky, S. Sharma, N. V. Dounskaia, "Directional biases reveal utilization of arm's biomechanical properties for optimization of motor behavior," Journal of Neurophysiology, vol. 98 no. 3, pp. 1240-1252, DOI: 10.1152/jn.00582.2007, 2007.

[11] K. Ohta, M. M. Svinin, Z. Luo, S. Hosoe, R. Laboissière, "Optimal trajectory formation of constrained human arm reaching movements," Biological Cybernetics, vol. 91 no. 1, pp. 23-36, 2004.

[12] Y. Zhang, J. Wang, Y. Xia, "A dual neural network for redundancy resolution of kinematically redundant manipulators subject to joint limits and joint velocity limits," IEEE Transactions on Neural Networks and Learning Systems, vol. 14 no. 3, pp. 658-667, DOI: 10.1109/TNN.2003.810607, 2003.

[13] L. Righetti, J. Buchli, M. Mistry, M. Kalakrishnan, S. Schaal, "Optimal distribution of contact forces with inverse-dynamics control," International Journal of Robotics Research, vol. 32 no. 3, pp. 280-298, DOI: 10.1177/0278364912469821, 2013.

[14] B. Siciliano, O. Khatib, Springer Handbook of Robotics,DOI: 10.1007/978-3-540-30301-5, 2008.

[15] T. Szkodny, A. Legowski, "The decreasing of 3D position errors in the system of two cameras by means of the interpolation method," vol. 9246, pp. 285-298, .

[16] Z. Nie, H. Xu, G. Yao, L. Wu, "Study on interpolation algorithm of printing screen transferred manipulator," Advanced Graphic Communications, Packaging Technology and Materials, vol. 369, pp. 675-681, 2016.

[17] M. Y. Kokurin, "Stable gradient projection method for nonlinear conditionally well-posed inverse problems," Journal of Inverse and ILL-Posed Problems, vol. 24 no. 3, pp. 323-332, DOI: 10.1515/jiip-2015-0047, 2016.

[18] J. Liu, C. Yang, J. Fang, "Improved optimization algorithm of the gram measurement matrix based on gradient projection," Journal of Huazhong University of Science and Technology: Nature Science, vol. 44 no. 8, pp. 62-65, 2016.

[19] W. Jin, Y. Censor, M. Jiang, "Bounded perturbation resilience of projected scaled gradient methods," Computational Optimization And Applications, vol. 63 no. 2, pp. 365-392, DOI: 10.1007/s10589-015-9777-x, 2016.

[20] B. Damas, J. Santos-Victor, "An online algorithm for simultaneously learning forward and inverse kinematics," Proceedings of the 25th IEEE/RSJ International Conference on Robotics and Intelligent Systems, IROS 2012, pp. 1499-1506, DOI: 10.1109/IROS.2012.6386156, .

[21] O. Kanoun, F. Lamiraux, P.-B. Wieber, "Kinematic control of redundant manipulators: Generalizing the task-priority framework to inequality task," IEEE Transactions on Robotics, vol. 27 no. 4, pp. 785-792, DOI: 10.1109/TRO.2011.2142450, 2011.

[22] C. Chen, D. Jackson, "Parameterization and evaluation of robotic orientation workspace: A geometric treatment," IEEE Transactions on Robotics, vol. 27 no. 4, pp. 656-663, DOI: 10.1109/TRO.2011.2135230, 2011.

[23] C. A. Pena, M. A. Guzman, P. F. Cardenas, "Inverse kinematics of a 6 dof industrial robot manipulator based on bio-inspired multi-objective optimization techniques," Proceedings of the IEEE Colombian Conference on Robotics and Automation (CCRA),DOI: 10.1109/CCRA.2016.7811428, .

[24] Z.-G. Hou, L. Cheng, M. Tan, "Multicriteria optimization for coordination of redundant robots using a dual neural network," IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics, vol. 40 no. 4, pp. 1075-1087, DOI: 10.1109/TSMCB.2009.2034073, 2010.

[25] Y. N. Zhang, "Analysis and design of recurrent neural networks and their applications to control and robotic systems," . 2003

[26] S. Ma, "A balancing technique to stabilize local torque optimization solution of redundant manipulators," Journal of Robotic Systems, vol. 13 no. 3, pp. 177-185, DOI: 10.1002/(SICI)1097-4563(199603)13:3<177::AID-ROB5>3.3.CO;2-I, 1996.

[27] Y. Zhang, D. Guo, S. Ma, "Different-level simultaneous minimization of joint-velocity and joint-torque for redundant robot manipulators," Journal of Intelligent & Robotic Systems, vol. 72 no. 3-4, pp. 301-323, 2013.

[28] Y. Zhang, X. Lv, Z. Li, Z. Yang, K. Chen, "Repetitive motion planning of PA10 robot arm subject to joint physical limits and using LVI-based primal-dual neural network," Mechatronics, vol. 18 no. 9, pp. 475-485, DOI: 10.1016/j.mechatronics.2008.04.005, 2008.

[29] Y. Zhang, W. Ma, X.-D. Li, "MATLAB simulink modeling and simulation of LVI-based primal-dual neural network for solving linear and quadratic programs," Neurocomputing, vol. 72 no. 7-9, pp. 1679-1687, DOI: 10.1016/j.neucom.2008.07.008, 2009.

[30] Y. Zhang, Z. Zhang, Repetitive Motion Planning and Control of Redundant Robot Manipulators, 2013.

[31] T. Zhang, H. An, H. Ma, "Joint torque and velocity optimization for a redundant leg of quadruped robot," International Journal of Advanced Robotic Systems, vol. 14 no. 5, 2017.

[32] M. J. Pérez-Ilzarbe, "Convergence analysis of a discrete-time recurrent neural network to perform quadratic real optimization with bound constraints," IEEE Transactions on Neural Networks and Learning Systems, vol. 9 no. 6, pp. 1344-1351, DOI: 10.1109/72.728385, 1998.

[33] M. J. Pérez-Ilzarbe, "New discrete-time recurrent neural network proposal for quadratic optimization with general linear constraints," IEEE Transactions on Neural Networks and Learning Systems, vol. 24 no. 2, pp. 322-328, DOI: 10.1109/TNNLS.2012.2223484, 2013.

[34] H. Tang, H. Li, Z. Yi, "A discrete-time neural network for optimization problems with hybrid constraints," IEEE Transactions on Neural Networks and Learning Systems, vol. 21 no. 7, pp. 1184-1189, DOI: 10.1109/TNN.2010.2049368, 2010.

[35] Q. Liu, J. Cao, "Global exponential stability of discrete-time recurrent neural network for solving quadratic programming problems subject to linear constraints," Neurocomputing, vol. 74 no. 17, pp. 3494-3501, DOI: 10.1016/j.neucom.2011.06.003, 2011.

[36] Y. Zhang, W. Ma, B. Cai, "From Zhang neural network to Newton iteration for matrix inversion," IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 56 no. 7, pp. 1405-1415, DOI: 10.1109/TCSI.2008.2007065, 2009.

[37] Y. Zhang, L. Xiao, G. Ruan, Z. Li, "Continuous and discrete time Zhang dynamics for time-varying 4th root finding," Numerical Algorithms, vol. 57 no. 1, pp. 35-51, DOI: 10.1007/s11075-010-9410-0, 2011.

[38] D. Guo, Y. Zhang, "Zhang neural network, Getz-Marsden dynamic system, and discrete-time algorithms for time-varying matrix inversion with application to robots' kinematic control," Neurocomputing, vol. 97, pp. 22-32, DOI: 10.1016/j.neucom.2012.05.012, 2012.

[39] Y. N. Zhang, Y. W. Yang, B. H. Cai, D. S. Guo, "Zhang neural network and its application to Newton iteration for matrix square root estimation," Neural Computing and Applications, vol. 21 no. 3, pp. 453-460, DOI: 10.1007/s00521-010-0445-x, 2012.

[40] Y. Zhang, B. Mu, H. Zheng, "Link between and comparison and combination of Zhang neural network and quasi-Newton BFGS method for time-varying quadratic minimization," IEEE Transactions on Cybernetics, vol. 43 no. 2, pp. 490-503, DOI: 10.1109/TSMCB.2012.2210038, 2013.

[41] L. Jin, Y. Zhang, "Discrete-time Zhang neural network of O( τ 3 ) pattern for time-varying matrix pseudoinversion with application to manipulator motion generation," Elsevier Science Publishers, 2014.

[42] B. Liao, Y. Zhang, L. Jin, "Taylor O(h 3 ) discretization of ZNN models for dynamic equality-constrained quadratic programming with application to manipulators," IEEE Transactions on Neural Networks and Learning Systems, vol. 27 no. 2, pp. 225-237, DOI: 10.1109/TNNLS.2015.2435014, 2016.

[43] L. Jin, Y. Zhang, "Discrete-time Zhang neural network for online time-varying nonlinear optimization with application to manipulator motion generation," IEEE Transactions on Neural Networks and Learning Systems, vol. 26 no. 7, pp. 1525-1531, DOI: 10.1109/TNNLS.2014.2342260, 2015.