2.1. Establishment of the workspace frame

As illustrated in Figure 1, a robotic additive manufacturing system is proposed, consisting of an industrial robot equipped with a printing end-effector, a rotating platform, and a binocular vision instrument (BVI). The robot is installed upside down on the gantry for additive manufacturing. Within the BVI coordinate system, the camera frame $\Sigma c$ is defined with its origin located at the midpoint between the two cameras and the z-axis aligned with the optical axis of the cameras. The coordinate frames of the robot base ( $\Sigma b$ ), robot wrist ( $\Sigma w$ ), and turning platform ( $\Sigma W$ ) are also established within the BVI software environment. The tool frame $\Sigma t$ is defined such that the x-axis is aligned with the direction connecting the two nozzles, and the z-axis is aligned with the wrist frame; the y-axis is determined by the right-hand rule. The relationship between the tool frame $\Sigma t$ on the end-effector and the wrist frame $\Sigma w$ can be obtained, as shown in Figure 2. The reflectors are mounted on the robot base, end-effector, and turntable, respectively, and are tracked by the BVI system. As a position-based visual feedback device, the BVI provides real-time pose information of the tracked objects using three or more non-collinear feature points. It is important to note that during any movement, all corresponding frames must remain within the field of view of the BVI.

Figure 1.

Robotic additive manufacturing system and space transformation of each frame.

Figure 2.

Schematic diagram of the tool frame.

Using VXelements software and the BVI (C-track 780) from Creamform, the positional and rotational information of the tracked model relative to the camera frame $\Sigma c$ can be acquired, recorded, and displayed simultaneously. Therefore, the homogeneous transform matrix ${}_b^c{\boldsymbol{T}}$ from $\Sigma b$ to $\Sigma c$ , ${}_w^c{\boldsymbol{T}}$ from $\Sigma w$ to $\Sigma c$ , and ${}_W^c{\boldsymbol{T}}$ from $\Sigma W$ to $\Sigma c$ can be obtained at any time. The homogeneous transformation matrix ${}_w^b{\boldsymbol{T}}$ , ${}_W^b{\boldsymbol{T}}$ and ${}_w^W{\boldsymbol{T}}$ can be obtained as follows

(1) \begin{align} {}_w^b{\boldsymbol{T}} &= {}_c^b{\boldsymbol{T}}{}_w^c{\boldsymbol{T}} = {({}_b^c{\boldsymbol{T}})^{ - 1}}{}_w^c{\boldsymbol{T}} \nonumber\\[4pt]{}_W^b{\boldsymbol{T}} &= {}_c^b{\boldsymbol{T}}{}_W^c{\boldsymbol{T}} = {({}_b^c{\boldsymbol{T}})^{ - 1}}{}_W^c{\boldsymbol{T}}\\[4pt]{}_w^W{\boldsymbol{T}} & = {}_c^W{\boldsymbol{T}}{}_w^c{\boldsymbol{T}} = {({}_W^c{\boldsymbol{T}})^{ - 1}}{}_w^c{\boldsymbol{T}}\nonumber \end{align}

Suppose $n$ feature points on the rigid end-effector are fixed and known from the definition of the tool frame $\Sigma w$ , whose homogeneous coordinates $ {}^w \boldsymbol{P}_i = \left[ \begin{array}{l@{\quad}l@{\quad}l@{\quad}l} {}^w{x_i}&{}^w{y_i}&{}^w{z_i}&1 \end{array} \right ]$ are denoted. With the pose $({x_t},{y_t},{z_t},{\alpha _t},{\beta _t},{\gamma _t})$ of $\Sigma w$ in $\Sigma c$ , the homogeneous transformation matrix ${}_w^c{\boldsymbol{T}}$ can be formulated as follows

(2) \begin{equation} { {}_w^c{\boldsymbol{T}} = \left [ {\begin{array}{c@{\quad}c} {R({\alpha _t},{\beta _t},{\gamma _t})}&{{{({x_t},{y_t},{z_t})}^T}}\\ 0&1 \end{array}} \right ]=\left [\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} c \alpha _t c \beta _t & c \alpha _t s \beta _t s \gamma _t-s \alpha _t c \gamma _t & c \alpha _t s \beta _t c \gamma _t+s \alpha _t s \gamma _t & x_t \\ s \alpha _t c \beta _t & s \alpha _t s \beta _t s \gamma _t+c \alpha _t c \gamma _t & s \alpha _t s \beta _t c \gamma _t-c \alpha _t s \gamma _t & y_t \\ -s \beta _t & c \beta _t s \gamma _t & c \beta _t c \gamma _t & z_t \\ 0 & 0 & 0 & 1 \end{array}\right ] } \end{equation}

where $R({\alpha _t},{\beta _t},{\gamma _t})$ is the rotation matrix from $\Sigma w$ in $\Sigma c$ and is taken from the ZYX-form Euler angle, $s$ and $c$ are the abbreviated forms of the functions $sin$ and $cos$ , respectively. Therefore, for the $i$ th feature point, the transformation equation as

(3) \begin{equation} {}^c{{\boldsymbol{P}}_i} = {}_w^c{\boldsymbol{T}}{}^w{{\boldsymbol{P}}_i} \end{equation}

The pose ${}^c{{\boldsymbol{P}}_i}$ of the manipulator can be obtained, with computations carried out in VXelements.

2.2. Trajectory planning algorithm for continuous path

The printing path is determined by the geometric shape of the additive manufacturing part, which consists of contour and infill paths. The slicing software divides the model into layers and generates cross-sectional paths based on the contours. For the printing plane with a spatial or inclined angle, the slicing software generates the infill path by setting the angle between the slicing plane and the printing base plane. A local interpolation frame is established in the slicing plane, and the coordinate transformation relationship between it and the printing base plane is obtained. In the local interpolation frame, the printing path is discretized into a series of planar free-form paths, which consist of straight lines and arcs. Therefore, trajectory interpolation for these straight lines and arcs is a critical aspect of robotic motion. In this paper, the trapezoidal velocity curve interpolation algorithm is employed for linear trajectory planning. This algorithm ensures smooth robotic movement, avoiding sudden speed changes that could cause shaking. For circular arc interpolation, a constant velocity interpolation method is used. Ensuring a smooth transition in velocity between the end of the linear segment and the circular arc. A schematic of linear and circular interpolation is shown in Figure 3.

Figure 3.

Diagram of linear and circular interpolation.

The NC code contains essential information such as robot pose, speed and acceleration limits, nozzle selection, and filament feed rate. When moving in continuous linear no-arc transition mode (the stop-at-turn motion), the robot is required to decelerate to a complete stop before changing direction and then re-accelerate afterward, which significantly increases the printing time [Reference Zhao, Zhu and Ding22]. Meanwhile, due to the filament material overflow caused by the nozzle pressure, the quality of the additive manufacturing products is further compromised when the robot stops at the corner. To enhance the speed continuity and print quality, the trajectory planning method for continuous linear arc transition is proposed. By limiting the maximum normal acceleration of small arcs, the speed at the corners can be restricted, and the accuracy effect is considered. The method results in smaller shape deviations and enables faster printing compared to the continuous linear no-arc transition method, thereby achieving a balance between print quality and manufacturing time.

For continuous short straight paths, a register is used to store the path information from the NC code for look-ahead planning, as illustrated in Figure 4. The register is first initialized, and the point data is sequentially stored from the $n-th$ container to the first container. The top of the register is then accessed to enable continuous linear arc transition planning, with the container data moving forward sequentially. After completing a segment of continuous path motion, a request is sent to the host computer system to retrieve a new set of data, thereby facilitating the update to the register.

Figure 4.

The design for the register.

Figure 5.

Continuous linear corners with arc transition.

For trajectory planning around corners, the robot’s velocity must transition smoothly at corner points to minimize speed fluctuations and improve motion stability. As discussed in Ref. [Reference Liu, Zhao, Li, Liu, Wang and Yue26], adjusting the robot’s speed to maintain continuity and smoothness through arc corners minimizes overflow and improves part quality. As shown in Figure 5, the additive manufacturing path employs smooth transitions using arcs, where the distance between the arc and the corner point is defined as the contour error $e_{co}$ . The deviation between the actual and the theoretical chord is termed the chord error $e_{ch}$ . Both errors are kept within acceptable limits. At the initial transition radius of the path at the $i-th$ corner point and the $(i+1)-th$ corner point, the radius is defined as follows

(4) \begin{equation} \begin{aligned} & r_i=\sin (\theta _i / 2) \cdot (e_{c o}-e_{c h}) / (1-\sin (\theta _i / 2)) \\ & r_{i+1}=\sin (\theta _{i+1} / 2) \cdot (e_{c o}-e_{c h}) / (1-\sin (\theta _{i+1} / 2)) \end{aligned} \end{equation}

where $\theta _i$ represents the angle between the two lines at the $i-th$ corner point. The projection lengths of the two arc segments on $\overrightarrow {{S_{i}}{S_{i+1}}}$ are

(5) \begin{equation} \begin{aligned} & L_{S_i P_{i,2}}=r_i / \tan (\theta _i / 2) \\ & L_{S_{i+1} P_{(i+1),1}}=r_{i+1} / \tan (\theta _{i+1} / 2) \end{aligned} \end{equation}

If $L_{S_i S_{i+1}}\lt L_{S_i P_{i,2}}+L_{S_{i+1} P_{(i+1),1}}$ the midpoint constraint method is required for initial adjustment.

(6) \begin{equation} \left \{\begin{array}{c} L_{S_i P_{i, 2}}=L_{S_i S_{i+1}}-L_{S_{i+1} P_{(i+1), 1}}, \min \left \{L_{S_i P_{i, 2}}, L_{S_{i+1} P_{(i+1), 1}}\right \}\lt 0.5 \times L_{S_i S_{i+1}} \cap L_{S_i P_{i, 2}}\gt L_{S_{i+1} P_{(i+1), 1}} \\ L_{S_{i+1} P_{i(i-1), 1}}=L_{S_i S_{i+1}}-L_{S_i P_{i, 2}}, \min \left \{L_{S_i P_{i, 2}}, L_{S_{i-1} P_{i(i+1), 1}}\right \}\lt 0.5 \times L_{S_i S_{i+1}} \cap L_{S_{i+1} P_{(i+1), 1}}\gt L_{S_i P_{i, 2}} \\ L_{S_i P_{i, 2}}=L_{S_{i+1} P_{i(i+1), 1}}=0.5 \times L_{S_i S_{i+1}}, \min \left \{L_{S_i P_{i, 2}}, L_{S_{i+1} P_{(i+1) 1}}\right \} \geq 0.5 \times L_{S_i S_{i+1}} \end{array}\right . \end{equation}

The corrected radius of the arc is given by

(7) \begin{equation} \begin{aligned} & r_i=L_{S_i P_{i,2}} \cdot \tan (\theta _i / 2) \\ & r_{i+1}=L_{S_{i+1} P_{(i+1), 1}} \cdot \tan (\theta _{i+1} / 2) \end{aligned} \end{equation}

The kinematic constraint for the motion along the arc segment is $v_{i c}^2=a \times r_i$ . However, the actual arc velocity is determined by the minimum value among the maximum allowable velocity, the maximum velocity under the chord error constraint, and the system’s rated maximum velocity, which is given by

(8) \begin{equation} v_c=\min (v_{i c}, v_{e_{-} c h}, v) \end{equation}

where $v_{e_{-} c h}=2 / t \times \sqrt {2 r \times e_{c h}-e_{c h}^2}$ , $t$ is interpolation period.

Furthermore, it is necessary to determine whether the linear trajectory can meet the velocity conditions. Specifically, this involves verifying whether the endpoint velocity of the previous arc can accelerate or decelerate to the starting velocity of the next arc. If this is not feasible, adjustments are made by fixing the arc radius and corner velocity on one side while modifying the arc radius and velocity on the other side, in conjunction with the velocity planning of the straight-line segment.

The detailed trajectory planning process is described as follows. $S_i$ represent the feature points in the toolpath trajectory, which are obtained from the NC code. A sharp corner is defined by two successive linear segments $\overrightarrow {{S_{i-1}}{S_{i}}}$ and $\overrightarrow {{S_{i}}{S_{i+1}}}$ of the toolpath, $a$ and $v$ represent the acceleration or deceleration and maximum speed, respectively. The small square points $A$ , $B$ , $C$ , and $D$ represent the speed change points. The robot trajectory typically follows a four-step velocity profile:

• • 1. Acceleration phase to reach the maximum target cruising speed, along the segment line $S_iA$ ,

• • 2. Uniform speed phase at constant speed $v_{i}$ , along the segment line $AB$ ,

• • 3. Deceleration phase while approaching the arc corner, along the segment line $BP_{i,1}$ ,

• • 4. Uniform speed phase, along the curve $P_{i,1}P_{i,2}$ .

Next, repeating the aforementioned steps. It is important to emphasize that if a line segment is short enough, it may bypass step 2, undergoing acceleration and direct deceleration. For example, the point $E$ on the right of Figure 5 serves as the transition point between acceleration and direct deceleration.

In the linear interpolation, let $L_{\mathrm{se}}=\left |P_{i, 2} P_{i+1,1}\right |$ be the length of the line segment, $v_s$ be the starting velocity, and $v_e$ be the ending velocity of the line segment, $v_{max}=v_i$ . If $L_{\mathrm{se}}\lt \left |\left (v_s^2-v_{\mathrm{e}}^2\right )\right | /(2a)$ , and $L_{\mathrm{se}}\lt \left |\left (v_s^2-v_{\mathrm{e}}^2\right )\right | /(2 a)$ , then there is no constant velocity stage in the trajectory, corresponding to Figure 6(b). The maximum velocity that the trajectory can reach is

(9) \begin{equation} v_{\lim }=\sqrt {L_{\mathrm{se}} \cdot a+\frac {v_{\mathrm{s}}^2+v_{\mathrm{e}}^2}{2}} \end{equation}

The acceleration duration $t_1$ and deceleration duration $t_2$ are respectively

(10) \begin{equation} t_1=\frac {v_{\max }-v_{\mathrm{s}}}{a}, t_2=\frac {v_{\max }-v_{\mathrm{e}}}{a} \end{equation}

The displacement at the interpolation time $t$

(11) \begin{equation} P_l(t)=\left \{\begin{array}{c@{\quad}c} v_s \cdot t+\dfrac {1}{2} a \cdot t^2, & 0 \leq t \leq t_1 \\[5pt] L_1+v_{\lim } \cdot t-\dfrac {1}{2} a \cdot t^2, & t_1\lt t \leq t_1+t_2 \end{array}\right . \end{equation}

where $L_1=v_s \cdot t_1+\left (a \cdot t_1^2\right ) / 2$ .

Otherwise, the trajectory includes acceleration, uniform speed, and deceleration stages, corresponding to Figure 6(a). In this case, the acceleration duration $t_1$ and deceleration duration $t_2$ are respectively

(12) \begin{equation} t_1=\frac {v_{\max }-v_{\mathrm{s}}}{a}, t_2=\frac {v_{\max }-v_{\mathrm{e}}}{a} \end{equation}

The uniform speed duration is

(13) \begin{equation} t_{\mathrm{m}}=\frac {L_{\mathrm{se}}-\left (\frac {v_{\max }^2-v_{\mathrm{s}}^2}{2 a}+\frac {v_{\max }^2-v_{\mathrm{e}}^2}{2 a}\right )}{v_{\max }} \end{equation}

The displacement at the interpolation time $t$

(14) \begin{equation} P_l(t)=\left \{\begin{array}{c@{\quad}c} v_s \cdot t+\dfrac {1}{2} a \cdot t^2, & 0 \leq t\lt t_1 \\ L_1+v_{\max } \cdot \left (t-t_1\right ), & t_1 \leq t\lt t_1+t_{\mathrm{m}} \\ L_2+v_{\max } \cdot \left (t-t_2\right )-\dfrac {1}{2} a \cdot \left (t-t_2\right ), & t_1+t_{\mathrm{m}} \leq t \leq t_1+t_2+t_{\mathrm{m}} \end{array}\right . \end{equation}

where $L_2=L_1+v_{\max } \cdot t_{\mathrm{m}}$ .

For circular arc interpolation, a uniform speed planning method is employed to derive the expressions of arc length and central angle as functions of time

(15) \begin{equation} \begin{aligned} &l(t)=v_e \cdot t\\ &\theta (t)=\theta \cdot l(t) / S \end{aligned} \end{equation}

where $S$ represents the arc length, $\theta$ denotes the central angle corresponding to the arc length. The interpolation points in the local coordinate system are expressed as a function of time

(16) \begin{equation} P_c(t)=[R \cdot \cos (\theta (t)), R \cdot \sin (\theta (t))] \end{equation}

Finally, the interpolation points are transformed from the local coordinate system to the printing base plane by the transformation matrix $_b^l T$ , yielding the interpolation points in the printing base plane is

(17) \begin{equation} \left \{\begin{array}{l} P(t)={ }_b^l T \cdot P_l(t) \text{ for line } \\ P(t)={ }_b^l T \cdot P_c(t) \text{ for arc } \end{array}\right . \end{equation}

Figure 6.

Schematic diagram of the phase trajectory: (a) Long-distance positioning; (b) Short-distance positioning.

Setting the velocity of the filament feeding motor to match the velocity of the robot’s end-effector, the velocity of the filament feeding motor $v_{fm}$ is

(18) \begin{equation} v_{f m}=\frac {v_a}{v} \cdot v_m \end{equation}

where $v$ is the maximum velocity for the robot in the NC code, $v_a$ is the actual speed of the robot’s end-effector, and $v_m$ is the rated speed of the feeding motor in the robot’s maximal motion velocity.

The execution logic of the planning algorithm is illustrated in Figure 7, encompassing NC code parsing, register storage, trajectory planning, and cross-interpolation of linear and circular trajectories. The sequential execution is indicated by the black arrow line, while parallel execution is represented by the red and green arrows. Notably, the red line signifies the arc interpolation input assignment, and the subsequent trajectory planning is executed in parallel with linear interpolation. The green line denotes the linear interpolation input assignment and the arc interpolation in parallel. Linear and arc interpolation are executed alternately. The pseudo-code for continuous linear arc transition is provided in Algorithm 1.

Figure 7.

Control flow chart for execution logic.

Algorithm 1.

Continuous linear arc transition mode based on the register, the planner, the interpolator, and then drive the robot.

In this algorithm, “M = 155” represents the continuous linear arc transition mode, $\boldsymbol{C}$ represents the register container (capacity of 10), and “planner (input)” and “planner (output)” represent inputs and outputs of the trajectory planner module, respectively. Similarly, the “ikine,” “line interpolator,” and “arc interpolator” represent the inverse kinematic solver, linear interpolator module, and circular interpolator module, respectively. The “para” represents an array of capacity 17 used to temporarily store the output of the planner. The “para [1:7]” represents the motion parameters of a line segment (start/end point coordinates and velocity), while the “para [8:17]” represents the motion parameters of an arc segment (start/mid/end point coordinates and velocity). $t_{line}$ and $t_{arc}$ represent the time of linear and circular interpolation.

2.3. Tracking control based on visual servoing

In the robotic operational space, the Cartesian-space vector $\boldsymbol{x}$ is defined as the position and orientation of the end-effector. According to the robot kinematic, the vector $\boldsymbol{x}$ is dependent on the joint configuration and can be expressed as a function of joint-space displacements as

(19) \begin{equation} {\boldsymbol{x}} = {\boldsymbol{h}}({\boldsymbol{q}}) \end{equation}

where ${\boldsymbol{h}}({\bullet})$ represents a nonlinear function that computes the operational space position variables from the joint space position variables. $\dot{\boldsymbol{x}}$ denotes a velocity vector in Cartesian space, which is related to the joint-space velocity vector $\dot{\boldsymbol{q}}$ by

(20) \begin{equation} {\dot{\boldsymbol{x}}} = {{\boldsymbol{J}}_m}({\boldsymbol{q}}){\dot{\boldsymbol{q}}} \end{equation}

where ${{\boldsymbol{J}}_m}({\boldsymbol{q}})$ is the link Jacobian matrix.

Consider a bounded continuous trajectory, which is designed far away enough from any singularity. In the joint space, the dynamics of a serial n-link rigid robot can be written as

(21) \begin{equation} \boldsymbol{M}(\boldsymbol{q})\ddot{\boldsymbol{q}} + \boldsymbol{C}(\boldsymbol{q,}\, \dot{\boldsymbol{q}})\dot{\boldsymbol{q}} + \boldsymbol{G}(\boldsymbol{q}) = \boldsymbol{\tau } \end{equation}

where $\boldsymbol{M(q)} \in \mathbb{R} {^{6 \times 6}}$ denotes the positive definite inertia matrix, $\boldsymbol{C}\boldsymbol{(}\boldsymbol{q,}\,\dot{\boldsymbol{q}}\boldsymbol{)} \in \mathbb{R} {^{6 \times 6}}$ is the centrifugal force and Coriolis force term, $\boldsymbol{G(q)} \in \mathbb{R} {^{6 \times 1}}$ denotes the gravity term matrix, and $\boldsymbol{\tau }$ is the joint control torque matrix.

The dynamic system in Cartesian space is the following form

(22) \begin{equation} {\boldsymbol{M}_x}(\boldsymbol{q})\ddot{\boldsymbol{x}} + {\boldsymbol{C}_x}(\boldsymbol{q,}\, \dot{\boldsymbol{q}})\dot{\boldsymbol{x}} + {\boldsymbol{G}_x}(\boldsymbol{q}) = {\boldsymbol{F}_x} \end{equation}

where ${{\boldsymbol{M}}_x}({\boldsymbol{q}}) = {\boldsymbol{J}}_m^{ - T}({\boldsymbol{q}}){\boldsymbol{M}}({\boldsymbol{q}}){\boldsymbol{J}}_m^{ - 1}({\boldsymbol{q}})$ , ${{\boldsymbol{C}}_x}({\boldsymbol{q}}) = {\boldsymbol{J}}_m^{ - T}({\boldsymbol{q}}){\boldsymbol{C}}({\boldsymbol{q}}){\boldsymbol{J}}_m^{ - 1}({\boldsymbol{q}})$ , ${{\boldsymbol{G}}_x}({\boldsymbol{q}}) = {\boldsymbol{J}}_m^{ - T}({\boldsymbol{q}}){\boldsymbol{G}}({\boldsymbol{q}})$ , ${{\boldsymbol{F}}_x} = {\boldsymbol{J}}_m^{ - T}({\boldsymbol{q}})\boldsymbol{\tau }$ .

$\textit {Property 1}$ : The centripetal-Coriolis matrix ${{\boldsymbol{C}}_x}({\boldsymbol{q}},{\dot{\boldsymbol{q}}})$ possesses the following property: ${{\boldsymbol{x}}^T}(\dot{\boldsymbol{M}}_x({\boldsymbol{q}}) - 2{{\boldsymbol{C}}_x}({\boldsymbol{q}},{\dot{\boldsymbol{q}}})){\boldsymbol{x}} = 0$ .

To meet the high precision requirements of additive manufacturing, it is necessary to compensate for the inherent inaccuracies of the robotic system. A visual servoing compensation strategy is established using the BVI, with measurements converted into pose information relative to the robot base frame. Feedback signals are returned to the servo compensation controller, which generates control commands to drive the robot, thereby adjusting its pose and achieving the desired position.

In the working space, $\boldsymbol{x}_d$ is the desired position in the desired trajectory, and $\dot {\boldsymbol{x}}_d$ and ${\ddot{\boldsymbol{x}}}_d$ represent the desired velocity and acceleration, respectively. Define the following equation

(23) \begin{equation} \begin{array}{l} \Delta {\boldsymbol{x}} = {{\boldsymbol{x}}_d} - {\boldsymbol{x}}\\ {{{\dot{\boldsymbol{x}}}}_r} = {{{\dot{\boldsymbol{x}}}}_d} - \textrm {T}\Delta {\boldsymbol{x}}\\ {\boldsymbol{r}} = {{{\dot{\boldsymbol{x}}}}_r} - {\dot{\boldsymbol{x}}} = \Delta {\dot{\boldsymbol{x}}} + \textrm { T}\Delta {\boldsymbol{x}} \end{array} \end{equation}

where, $\Delta {\boldsymbol{x}}$ is the pose error in Cartesian space, and $\boldsymbol{x}_d$ and $\dot {\boldsymbol{x}}_d$ are obtained by the interpolation module. In the Cartesian-space controller, the Cartesian-space feedback $\boldsymbol{x}$ and $\dot {\boldsymbol{x}}$ are usually measured by the BVI, which undergoes image processing, filtering, and coordinate transformation. However, due to the presence of measurement noise, directly differentiating position values $\boldsymbol{x}$ results in significant errors in velocity $\dot {\boldsymbol{x}}$ . To mitigate this issue, a velocity observer is introduced to estimate velocity using measured position data. $\mathrm{T}$ is a positive definite matrix.

The dynamic parameters were predicted using the RBF neural networks. The controller in Cartesian space is designed as

(24) \begin{equation} {{\boldsymbol{F}}_x} = \hat{\boldsymbol{w}}_M^T \cdot {{\boldsymbol{h}}_M}({\boldsymbol{q}}){{\ddot{\boldsymbol{x}}}_r} + \hat{\boldsymbol{w}}_C^T \cdot {{\boldsymbol{h}}_C}({\boldsymbol{z}}){{\dot{\boldsymbol{x}}}_r} + \hat{\boldsymbol{w}}_G^T \cdot {{\boldsymbol{h}}_G}({\boldsymbol{q}}) + {\boldsymbol{Kr}} + {k_s}\mathop {\textrm {sgn}} ({\boldsymbol{r}}) \end{equation}

where $(\hat{\bullet})$ is the estimation of $({\bullet})$ , the first three terms of the control law are model-based control components predicted by the RBF neural networks, and the construction process is detailed in Appendix A. The $\boldsymbol{Kr}$ is the proportional derivative (PD) control term, ${{\boldsymbol{K}}_p} = {\boldsymbol{K}}\textrm {T}$ and ${{\boldsymbol{K}}_d} = {\boldsymbol{K}}$ are diagonal and positive definite matrices, making ${\textrm {T}_{ij}} = \frac {1}{{{K_{ij}}^2}}(i = j)$ . The last term is a robust term to suppress the neural networks’ modeling errors, and ${k_s} \gt \left \| {\boldsymbol{E}} \right \|$ .

The control parameters $\boldsymbol{K}$ corresponding to the four stages in Figure 6 are designed as

(25) \begin{equation} \left \{ {\begin{array}{l} {{{\boldsymbol{K}}_{\boldsymbol{I}}} = {\bar{\boldsymbol{K}}} - \dfrac {1}{{2{a_1}}}\Delta {\boldsymbol{x}} \cdot {\boldsymbol{I}}}\\ {{{\boldsymbol{K}}_{\boldsymbol{II}}} = {\bar{\boldsymbol{K}}}}\\ {{{\boldsymbol{K}}_{\boldsymbol{III}}} = {\bar{\boldsymbol{K}}} + \dfrac {1}{{2{a_2}}}\Delta {\boldsymbol{x}} \cdot {\boldsymbol{I}}}\\ {{{\boldsymbol{K}}_{\boldsymbol{IV}}} = \dfrac {{{\bar{\boldsymbol{K}}}}}{2}} \end{array}} \right . \end{equation}

where $\bar{\boldsymbol{K}}$ is the nominal value, $a_1$ and $a_2$ represent the value of acceleration and deceleration, respectively. $\boldsymbol{I}$ represents the unit matrix. In the constant acceleration phase, parameters are designed to achieve a smooth and rapid start of the motor. Once the velocity reaches its maximum, the control parameters constrain the operating speed. During deceleration, the parameters ensure a stable and accurate approach to the target position. For short-distance positioning tasks, the uniform speed phase (step 2) is omitted.

The block diagram of the visual servoing controller integrated with a velocity observer is illustrated in Figure 8. The matrix ${{\boldsymbol{J}}_m}({\boldsymbol{q}})$ denotes the Jacobian matrix, and the feedback information in Cartesian space is directly transformed into the controller output of the robot in joint space using the Jacobian matrix. $\boldsymbol{\tau }$ is the torque exerted in the joint space.

Figure 8.

The principle of the visual servoing controller.

Define the adaptive law of neural networks weights [Reference Wei, Li, Tian, Wei, Bai, Zhang, Li and Hou29]:

(26) \begin{equation} \begin{array}{l} \dot{\hat{\boldsymbol{W}}}_{Mk} = {{\boldsymbol{\Gamma }}_{Mk}} \cdot {\boldsymbol{\phi } _{Mk}}({\boldsymbol{q}}){{{\ddot{\boldsymbol{x}}}}_r}{{\boldsymbol{r}}_k}\\ \dot{\hat{\boldsymbol{W}}}_{Ck} = {{\boldsymbol{\Gamma }}_{Ck}} \cdot {\boldsymbol{\phi } _{Ck}}({\boldsymbol{z}}){{{\dot{\boldsymbol{x}}}}_r}{{\boldsymbol{r}}_k}\\ \dot{\hat{\boldsymbol{W}}}_{Gk} = {{\boldsymbol{\Gamma }}_{Gk}} \cdot {\boldsymbol{\phi } _{Gk}}({\boldsymbol{q}}){{\boldsymbol{r}}_k} \end{array} \end{equation}

where ${\boldsymbol{\Gamma }}_{MK}$ , ${\boldsymbol{\Gamma }}_{CK}$ , and ${\boldsymbol{\Gamma }}_{GK}$ are the symmetric positive definite matrix, and ${\boldsymbol{\phi } _{Mk}}({\boldsymbol{q}})$ , ${\boldsymbol{\phi } _{Ck}}({\boldsymbol{z}})$ , and ${\boldsymbol{\phi } _{Gk}}({\boldsymbol{q}})$ are the Gaussian kernel functions of the neural networks, respectively.

The velocity observer is proposed to solve the velocity estimation

(27) \begin{equation} { \begin{array}{l} {\dot {\hat {\boldsymbol{x}}}} = \boldsymbol{P} + \Lambda {\tilde {\boldsymbol{x}}}\\ \dot{\boldsymbol{p}} = \boldsymbol{M}_x^{ - 1}(\boldsymbol{q})[{\boldsymbol{F}_x} - {\boldsymbol{C}_x}(\boldsymbol{q,}\, \dot{\boldsymbol{q}}) \dot {\hat {\boldsymbol{x}}} - {\boldsymbol{G}_x}(\boldsymbol{q})] \end{array} } \end{equation}

where $\boldsymbol{p}$ is an auxiliary variable, $\Lambda$ is a diagonal positive definite matrix designed based on the stability requirements of the closed-loop system. The expression $\Lambda$ will be presented in the stability analysis of the observer system. The dynamic coefficients are parameter values predicted using the RBF neural networks. Note that the equation is basically the robot Eq. (22) with $\dot {\boldsymbol{x}}$ replaced by $\dot {\hat {\boldsymbol{x}}}$ .

Before the stability of the observer systems is analyzed, the observation error $\dot {\tilde {\boldsymbol{x}}}$ for the system is derived. By taking the time derivative of $\dot {\hat {\boldsymbol{x}}}$

(28) \begin{equation} {{\ddot {\hat {\boldsymbol{x}}}}} = \dot {\boldsymbol{p}} + \Lambda {{\dot {\tilde {\boldsymbol{x}}}}} \end{equation}

After pre-multiplying Eq. (28) by ${\boldsymbol{M}_x}(\boldsymbol{q})$ , and substituting $\dot {\boldsymbol{p}}$ of Eq. (27)

(29) \begin{equation} {\boldsymbol{M}_x}(\boldsymbol{q}){\ddot {\hat {\boldsymbol{x}}}} = {\boldsymbol{F}_x} - {\boldsymbol{C}_x}(\boldsymbol{q,}\, \dot{\boldsymbol{q}})\dot {\hat {\boldsymbol{x}}} - {\boldsymbol{G}_x}(\boldsymbol{q}) + \Lambda {\boldsymbol{M}_x}(\boldsymbol{q})\dot {\tilde {\boldsymbol{x}}} \end{equation}

Subtracting Eq. (29) from Eq. (22), we obtain the closed-loop error system for Cartesian-space velocity observation error $\dot {\tilde {\boldsymbol{x}}}$ as

(30) \begin{equation} {\boldsymbol{M}_x}(\boldsymbol{q}){\ddot {\tilde {\boldsymbol{x}}}} = - {\boldsymbol{C}_x}(\boldsymbol{q,}\, \dot{\boldsymbol{q}})\dot {\tilde {\boldsymbol{x}}} - \Lambda {\boldsymbol{M}_x}(\boldsymbol{q})\dot {\tilde {\boldsymbol{x}}} \end{equation}

Let the function $V_0$ be defined as follows to analyze the velocity observer stability.

(31) \begin{equation} {V_0} = \frac {1}{2}{{\dot {\tilde {\boldsymbol{x}}}}^T}{{\boldsymbol{M}}_x}({\boldsymbol{q}}){\dot {\tilde {\boldsymbol{x}}}} \end{equation}

Taking the time derivative of $V_0$

(32) \begin{equation} {\dot {V}_0} = {{\dot {\tilde {\boldsymbol{x}}}}^T}{{\boldsymbol{M}}_x}({\boldsymbol{q}}){\ddot {\tilde {\boldsymbol{x}}}} + \frac {1}{2}{{\dot {\tilde {\boldsymbol{x}}}}^T}{{\dot {\boldsymbol{M}}}_x}({\boldsymbol{q}}){\dot {\tilde {\boldsymbol{x}}}} \end{equation}

After substituting the error system for $\dot {\tilde {\boldsymbol{x}}}$ of Eq. (30) from Eq. (32), results in

(33) \begin{equation} \begin{split} {\dot {V}_0} &= \frac {1}{2}{{\dot {\tilde {\boldsymbol{x}}}}^T}{{\dot {\boldsymbol{M}}}_x}({\boldsymbol{q}}){\dot {\tilde {\boldsymbol{x}}}} + {{{\dot {\tilde {\boldsymbol{x}}}}}^T}( - {{\boldsymbol{C}}_x}({\boldsymbol{q}},{\dot{\boldsymbol{q}}}){\dot{\tilde{\boldsymbol{x}}}} - \Lambda {{\boldsymbol{M}}_x}({\boldsymbol{q}}){\dot {\tilde {\boldsymbol{x}}}})\\ &=\frac {1}{2}{{\dot {\tilde {\boldsymbol{x}}}}^T}({{\dot {\boldsymbol{M}}}_x}({\boldsymbol{q}})-2{{\boldsymbol{C}}_x}({\boldsymbol{q}},{\dot{\boldsymbol{q}}})){\dot{\tilde{\boldsymbol{x}}}}-{\dot {\tilde {\boldsymbol{x}}}}^T \Lambda {{\boldsymbol{M}}_x}({\boldsymbol{q}}){\dot {\tilde {\boldsymbol{x}}}} \end{split} \end{equation}

Applying the $\textit {property 1}$ , we obtain

(34) \begin{equation} {\dot {V}_0}=-{\dot {\tilde {\boldsymbol{x}}}}^T \Lambda {{\boldsymbol{M}}_x}({\boldsymbol{q}}){\dot {\tilde {\boldsymbol{x}}}} \le - \Lambda {m_1} \| {\dot {\tilde {\boldsymbol{x}}}} \| \le 0 \end{equation}

$V_0$ is positive definite, $\dot {V_0}$ is semi-negative definite, and the closed-loop system has progressive stabilization. When $t \to \infty$ , ${\dot {\tilde {\boldsymbol{x}}}} \to \infty$ , the observer serves as an estimator of velocity effectively.

The robot dynamics Eq. (22) is written as Eq. (35) according to Appendix A.

(35) \begin{equation} \begin{array}{l} {{\boldsymbol{F}}_x} = \left ( {{{\boldsymbol{F}}_{MNN}}({\boldsymbol{q}}) + {{\boldsymbol{E}}_M}} \right ){\ddot{\boldsymbol{x}}} + \left ( {{{\boldsymbol{F}}_{CNN}}({\boldsymbol{z}}) + {{\boldsymbol{E}}_{\boldsymbol{C}}}} \right ){\dot{\boldsymbol{x}}} + \left ( {{{\boldsymbol{F}}_{GNN}}({\boldsymbol{q}}) + {{\boldsymbol{E}}_{\boldsymbol{G}}}} \right ) + {k_s}\mathop {\textrm {sgn}} ({\boldsymbol{r}})\\ = \left [ {{\boldsymbol{w}}_M^T{{\boldsymbol{h}}_M}({\boldsymbol{q}}) + {{\boldsymbol{E}}_M}} \right ]{\ddot{\boldsymbol{x}}} + \left [ {{\boldsymbol{w}}_C^T{{\boldsymbol{h}}_C}({\boldsymbol{z}}) + {{\boldsymbol{E}}_C}} \right ]{\dot{\boldsymbol{x}}} + \left [ {{\boldsymbol{w}}_G^T{{\boldsymbol{h}}_G}({\boldsymbol{q}}) + {{\boldsymbol{E}}_G}} \right ] + {k_s}\mathop {\textrm { sgn}} ({\boldsymbol{r}}) \end{array} \end{equation}

Applying the control law from Eq. (24) to the robot dynamics in Eq. (35) and using Eq.(23) yields the tracking error equation

(36) \begin{equation} \begin{aligned} & \left [ \boldsymbol{w}_M^T \boldsymbol{h}_M (\boldsymbol{q}) + \boldsymbol{E}_M \right ]{\ddot{\boldsymbol{r}}} + \left [ \boldsymbol{w}_C^T \boldsymbol{h}_C({\boldsymbol{z}}) + \boldsymbol{E}_C \right ]\dot{\boldsymbol{r}}\\ &= {\boldsymbol{M} \dot{\boldsymbol{r}}} + {\boldsymbol{Cr}}= {\tilde {\boldsymbol{w}}}_M^T{{\boldsymbol{h}}_M}({\boldsymbol{q}})\ddot {\boldsymbol{x}}_r + {\tilde {\boldsymbol{w}}}_C^T{{\boldsymbol{h}}_C}({\boldsymbol{z}})\dot {\boldsymbol{x}}_r + {\tilde {\boldsymbol{w}}}_G^T{{\boldsymbol{h}}_G}({\boldsymbol{q}}) + {\boldsymbol{E}} - {\boldsymbol{Kr}} - {k_s}\mathop {\textrm {sgn}} ({\boldsymbol{r}}) \end{aligned} \end{equation}

where $(\tilde{\bullet})=({\bullet})-(\hat{\bullet})$ represents the estimation error of the parameter, ${\boldsymbol{E}} = {{\boldsymbol{E}}_M}{{\ddot{\boldsymbol{x}}}_r} + {{\boldsymbol{E}}_C}{{\dot{\boldsymbol{x}}}_r} + {{\boldsymbol{E}}_G}$ .

To prove the stability of the system predicted by the RBF neural networks, the Lyapunov function is defined as

(37) \begin{equation} {V_1} = \frac {1}{2}{{\boldsymbol{r}}^T}{{\boldsymbol{M}}_x}{\boldsymbol{r}} + \frac {1}{2}\sum \limits _{k = 1}^n \tilde{\boldsymbol{w}}_{Mk}^T {\boldsymbol{\Gamma }}_{Mk}^{ - 1}\tilde{\boldsymbol{w}}_{Mk} + \frac {1}{2}\sum \limits _{k = 1}^n \tilde{\boldsymbol{w}}_{Ck}^T {\boldsymbol{\Gamma }}_{Ck}^{ - 1}\tilde{\boldsymbol{w}}_{Ck} + \frac {1}{2}\sum \limits _{k = 1}^n \tilde{\boldsymbol{w}}_{Gk}^T {\boldsymbol{\Gamma }}_{Gk}^{ - 1}\tilde{\boldsymbol{w}}_{Gk} \end{equation}

Computing the derivative of Eq. (37) and simplifying it to get

(38) \begin{align} \dot{V}_1 = {\boldsymbol{r}}^T \boldsymbol{M} \dot{\boldsymbol{r}} + \frac {1}{2}\boldsymbol{r}^T {\dot{\boldsymbol{Mr}}} + \sum \limits _{k = 1}^n \tilde{\boldsymbol{w}}_{Mk}^T {\boldsymbol{\Gamma }}_{Mk}^{ - 1}\dot{\tilde{\boldsymbol{w}}}_{Mk} + \sum \limits _{k = 1}^n \tilde{\boldsymbol{w}}_{Ck}^T {\boldsymbol{\Gamma }}_{Ck}^{ - 1}{\dot{\tilde{\boldsymbol{w}}}_{Ck}} + \sum \limits _{k = 1}^n \tilde{\boldsymbol{w}}_{Gk}^T {\boldsymbol{\Gamma }}_{Gk}^{ - 1}{\dot{\tilde{\boldsymbol{w}}}}_{Gk} \end{align}

According to $\textit {property 1}$ , and $ \dot{\tilde{\boldsymbol{w}}}_{MK} = - \dot{\hat{\boldsymbol{w}}}_{Mk}$ , $ \dot{\tilde{\boldsymbol{w}}}_{CK} = -{\dot{\hat{\boldsymbol{w}}}}_{Ck}$ , $ \dot{\tilde{\boldsymbol{w}}}_{GK} = -\dot{\hat{\boldsymbol{w}}}_{Gk}$ . Substituted into the above formula

(39) \begin{equation} {\dot V_1} = {{\boldsymbol{r}}^T}(\boldsymbol{M}\dot{\boldsymbol{r}} + {\boldsymbol{Cr}}) - \sum \limits _{k = 1}^n \tilde{\boldsymbol{w}}_{Mk}^T {\boldsymbol{\Gamma }}_{Mk}^{ - 1}{\dot{\hat{\boldsymbol{w}}}}_{Mk} - \sum \limits _{k = 1}^n \tilde{\boldsymbol{w}}_{Ck}^T {\boldsymbol{\Gamma }}_{Ck}^{ - 1}{\dot{\hat{\boldsymbol{w}}}}_{Ck} - \sum \limits _{k = 1}^n \tilde{\boldsymbol{w}}_{Gk}^T {\boldsymbol{\Gamma }}_{Gk}^{ - 1}{\dot{\hat{\boldsymbol{w}}}_{Gk}} \end{equation}

Taking Eq. (36) into Eq. (39)

(40) \begin{equation} \begin{array}{l} {{\dot V}_1} = - {{\boldsymbol{r}}^T}{\boldsymbol{Kr}} + {{\boldsymbol{r}}^T}{\boldsymbol{E}} - {k_s}{{\boldsymbol{r}}^T}\mathop {\textrm {sgn}} ({\boldsymbol{r}}) + {{\boldsymbol{r}}^T}\tilde{\boldsymbol{w}}_M^T \cdot {{\boldsymbol{h}}_M}({\boldsymbol{q}}){{{\ddot{\boldsymbol{x}}}}_r} + {{\boldsymbol{r}}^T}\tilde{\boldsymbol{w}}_C^T \cdot {{\boldsymbol{h}}_C}({\boldsymbol{z}}){{{\dot{\boldsymbol{x}}}}_r} + {{\boldsymbol{r}}^T}\tilde{\boldsymbol{w}}_G^T \cdot {{\boldsymbol{h}}_G}({\boldsymbol{q}})\\ - \sum \limits _{k = 1}^n {\tilde{\boldsymbol{w}}_{Mk}^T} {\boldsymbol{\Gamma }}_{Mk}^{ - 1}{{\dot{\hat{\boldsymbol{w}}}}_{Mk}} - \sum \limits _{k = 1}^n {\tilde{\boldsymbol{w}}_{Ck}^T} {\boldsymbol{\Gamma }}_{Ck}^{ - 1}{{\dot{\hat{\boldsymbol{w}}}}_{Ck}} - \sum \limits _{k = 1}^n {\tilde{\boldsymbol{w}}_{Gk}^T} {\boldsymbol{\Gamma }}_{Gk}^{ - 1}{{\dot{\hat{\boldsymbol{w}}}}_{Gk}} \end{array} \end{equation}

Substitute the adaptive rate Eq. (26) into the above equation and ${k_s} \gt \| {\boldsymbol{E}} \|$ , we obtain

(41) \begin{equation} {\dot V_1} = - {{\boldsymbol{r}}^T}{\boldsymbol{Kr}} + {{\boldsymbol{r}}^T}{\boldsymbol{E}} - {k_s}{{\boldsymbol{r}}^T}\mathop {\textrm {sgn}} ({\boldsymbol{r}}) \le 0 \end{equation}

In general, ${V_1} \ge 0$ , ${\dot {V_1}} \le 0$ , when the manipulator operates in a workspace where the kinematics mapping between Cartesian space and joint space is one-to-one, and the Jacobian matrix is nonsingular, the asymptotic stability can be concluded according to LaSalle’s Invariance Theorem.