1. Introduction
Cable-driven parallel robots (referred to as CDPRs or CDPR hereafter) are special parallel robots that use cables instead of rigid links to transmit motion and forces [Reference Mamidi and Bandyopadhyay1]. The flexible nature of cable-driven robots provides advantages over traditional parallel robots, such as low moving inertia, large workspace, high velocity, and low maintenance costs [Reference Qian, Zhao, Qian, Wang and Zi2]. As a result, they have attracted considerable attention from researchers and have been widely applied in astronaut training [Reference Zhang, Zou, Wang and Pei3], rehabilitation robotics [Reference Khadem, Inel, Carbone and Tich4], giant radio telescopes [Reference Duan, Qiu, Zhang and Zi5], and human-machine interaction devices [Reference Fortin-Côté, Cardou and Gosselin6].
The most basic model of CDPRs is the geometric model, which focuses on the robot’s geometric layout and kinematics while simplifying the dynamics. This model traditionally assumes that the cables are idealized as precise linear distances between two points in space, specifically between the proximal anchor points on the fixed frame and the distal anchor points on the moving platform, with the distances perfectly matched. A common extension to the standard model (integrating geometry, dynamics, and control) is the inclusion of pulleys and friction [Reference Kraus, Kessler and Pott7–Reference Liu, Li, Zhang, Hu and Zhang9]. For example, Mamidi et al. incorporated pulley kinematics based on the recursive sub-system-level Lagrangian multiplier approach (RSSLM) [Reference Mamidi and Bandyopadhyay10]. Fabritius et al. proposed a framework aiming to unify various CDPR models found in the literature. It consists of four nested CDPR model stages (point platform, geometric, elastic, and catenary) and two optional model extensions (force offset and pulleys) [Reference Fabritius, Miermeister, Kraus and Pott11]. Patel et al. described a CDPR end effector considering pulleys by using closed-loop equations for cable lengths and cable length squares, as well as an unconstrained parameterization of the pose [Reference Patel, Nguyen and Caverly12].
Due to the non-negligible mass and elasticity of the cables, the cables deform when the end effector moves [Reference Mamidi and Bandyopadhyay1, Reference Peng, Xiao, Chen, Wei, Lin and Zhuo13]. Therefore, in motion planning, it is more inclined to study the tension distribution of cables rather than their lengths. Since cables can only exert tension, redundancy is required to fully constrain the system [Reference Abbasnejad, Eden and Lau14, Reference Ding, Zheng, Liu, Guo and Guo15]. In other words, the number of cables
$m$
in CDPRs is greater than the number of degrees of freedom
$n$
of the end effector. This configuration is referred to as overconstrained, and the redundancy is denoted as
$r = m - n$
. Typically, the redundancy is set to 1–3 [Reference Ferravante, Riva, Taghavi, Braghin and Bock16]. However, for overconstrained systems, since
$m$
is not equal to
$n$
, the number of equations for cable tensions is smaller than the number of variables, resulting in infinitely many solutions for cable tensions [Reference Sun, Liu, Gao, Li, Ding and Deng17]. At time
$t$
, the set of all solutions forms an
$r$
-dimensional linear subspace. The problem of force distribution in redundantly constrained CDPRs is inevitable and often time-consuming in control. Chawla et al. trained a neural network for the inverse statics (IKS) problem of CDPRs, providing solutions for both inverse statics and forward statics and significantly improving computational speed [Reference Chawla, Pathak, Notash, Samantaray, Li and Sharma18]. Martin-Parra et al. proposed a new design scheme, established the kinematic and static models of a new robot, verified the static model of the robot through an optimized interior point algorithm, and obtained cable force distributions for different end effector trajectories [Reference Martin-Parra, Juarez-Perez, Gonzalez-Rodriguez, Gonzalez-Rodriguez, Lopez-Diaz and Rubio-Gomez19].
In order to obtain stable and continuous tension solutions (cable tension exhibits smooth temporal variation without abrupt changes) to address the redundancy problem in CDPRs, many scholars have constructed an objective function that treats the dynamics or statics equations of CDPRs as constraints, thereby transforming the problem of tension distribution of flexible cables into an optimization problem. Oh et al. utilized the 1-norm of the tension set as the objective function and provided methods for solving linear programming and quadratic programming problems to obtain feasible solutions within the operational space of the system [Reference Oh and Agrawal20]. Borgstrom proposed a new linear programming formulation that generates the optimal and safe tension distribution in CDPRs by introducing relaxation variables [Reference Borgstrom, Jordan, Sukhatme, Batalin and Kaiser21]. However, the solution obtained through 1-norm optimization is not continuous and may cause high-frequency vibrations [Reference Bruckmann, Pott and Hiller22]. Taghirad et al. introduced the 2-norm of the tension set as the objective function for solving the problem [Reference Taghirad and Bedoustani23], while Gosselin et al. introduced the p-norm as the optimization objective [Reference Gosselin and Grenier24]. Geng proposed an optimization algorithm based on the minimum tension-to-hyperstaticity ratio to address excessive stiffness [Reference Geng, Li, Liu, Li, Zheng and Li25]. However, most of these optimization algorithms are iterative, resulting in tension sets that often approach the lower limit of cable tension and lack real-time performance. Pott et al. proposed an iterative algorithm solely based on matrix operations [Reference Pott26, Reference Pott and Bruckmann27]. Fabritius et al. proposed a force correction method based on the null space to improve the dynamic performance of CDPRs. This method modifies the cable forces within the null space of the robot’s structural matrix, keeping them within feasible limits and as close as possible to the desired level without altering the platform’s posture [Reference Fabritius, Rubio-Gómez, Martin, Santos, Kraus and Pott28]. Jin et al. introduced a modified index to evaluate the motion/force transmission performance of CDPRs in order to address the tension constraint of flexible cables. By calculating the power coefficients of positive and negative output motion directions, they ensured that the driving forces in the opposite directions of the two motion screws are positive. Furthermore, through optimization design based on the transmission index, CDPRs can reduce singularities in the workspace [Reference Jin, Ye and Li29]. Einar Ueland et al. applied the logarithmic barrier method to construct the objective function and demonstrated its continuous differentiability with respect to the cable tensions at different time instances when it obtains the optimal solution [Reference Ueland, Sauder and Skjetne30]. They also introduced relaxation variables to enable rapid convergence of the proposed optimization problem, resulting in better real-time performance compared to other objective functions.
Due to the iterative nature of tension distribution algorithms based on optimization problems, geometric methods have attracted considerable attention from scholars. Marc Gouttefarde et al. introduced an algorithm for the tension feasible region (TFR) of n-DOF (degrees of freedom) CDPRs specifically designed for n + 2 cables driving n platforms. The algorithm determines the vertices of the 2-D convex polygon in either clockwise or counterclockwise direction and involves corresponding improvements to determine various cable tension distributions [Reference Gouttefarde, Lamaury, Reichert and Bruckmann31]. Cui et al. proposed a non-iterative optimization method for cable tension using the Graham scan method incorporated into the geometric approach. It was applied in a dual-thread PID force/position hybrid controller, providing a non-iterative method for optimizing the cable tension of CDPRs with
$r=2$
, simultaneously addressing computational efficiency, minimizing cable tension, optimizing safety, and stiffness concerns [Reference Cui, Tang, Hou and Sun32]. Rodriguez-Barroso et al. presented a centroid-based approach to compute the polyhedron of cable tension possibilities, thereby calculating the maximum torque that can be exerted by a suspended cable-driven robot and determining the optimal tension distribution required to achieve any desired maximum torque [Reference Rodriguez-Barroso and Saltaren33]. Sun et al. developed a geometric tension feasible region algorithm by interpreting the TFR as the convex intersection of multiple boundary-parallel spaces, computing the tension feasible region for CDPRs with redundancy
$r=2$
and 3 [Reference Sun, Liu, Gao, Li, Ding and Deng17]. Mohammad Reza Mousavi et al. proposed an algorithm applicable to redundancy levels
$r=1-3$
, aiming to keep the cables within their tensile limits, preventing tearing during operations, defining limits on cable forces using a problem formulation, calculating the optimal solution for cable tension, tracking paths, and ensuring the safety of all actuators [Reference Mousavi, Ghanbari, Moosavian and Zarafshan34]. Geometric methods have several advantages, such as short computation time, applicability within the entire feasible workspace of force screws, and wide usage. They are the most widely used and popular non-iterative methods. With a good initial value, the computational cost of geometric methods can be significantly reduced. However, logical judgments are employed in each iteration step of geometric methods. Therefore, geometric methods cannot provide an explicit function for the optimal set of tensions [Reference Geng, Li, Liu, Li, Zheng and Li25]. Meanwhile, Einar Ueland et al. mentioned that the geometric algorithm involves the calculation of the 2-norm and pseudoinverse when calculating cable tension, which may result in discontinuous tension values [Reference Ueland, Sauder and Skjetne30].
This paper proposes a new force distribution algorithm for a CDPR with four cables and 3-DOF. Although many optimized iterative algorithms have been proposed to solve the redundancy constraint force distribution problem for CDPRs, issues with poor real-time performance have not been fully resolved, and there are certain requirements for hardware conditions. Additionally, geometric methods, when the redundancy
$r$
is 1, only use the geometric boundary of line segments for tension distribution. Due to insufficient research, effective tension analysis methods have not yet been provided for CDPRs with a redundancy of
$r=1$
. This study begins with the wave equation from elastic dynamics to derive the coefficients that influence force distribution, enabling the distribution of cable forces without iterations and without the need for a complete calculation of the Jacobian matrix. Moreover, this paper investigates the impact of primary coefficients on cable tension, ensuring the continuity of tension. Lastly, a specific numerical example is presented to simulate and analyze the algorithm’s real-time performance, the continuity, and the accuracy of the cable tension.
The organization of the paper is as follows: Section 2 introduces the wave equation. Section 3 starts from the general solution of the wave equation and derives an integral equation about the tension distribution coefficient of the cable, and proves the existence of a solution in the form of a cubic polynomial. Section 4 gives a specific numerical example to verify the correctness of the proposed algorithm, demonstrates its real-time performance, and compares it with existing methods. Section 5 summarizes the main features of this algorithm.
2. Introduction to the wave equation
The wave equation is a hyperbolic partial differential equation that describes wave propagation in continuum mechanics and physics. For an
$n$
-dimensional continuous medium, under zero initial conditions, the wave equation can be expressed as follows:
\begin{equation} \begin{split} & \Box ^2 \phi = -g(x) \\ & \phi (\boldsymbol{x},0) = \dot{\phi } (\boldsymbol{x},0) = 0 \end{split} \end{equation}
Wherein,
$\boldsymbol{x}=\left [ x_1,x_2,\cdots, x_n \right ]$
represents the various axes/dimensions of the object of study. Its general solution is given by:
\begin{equation} \begin{split} \phi (\boldsymbol{x}, t) &= \\ &\frac{c_d}{2\pi \rho } \int _C \frac{g(\boldsymbol{x}', t')}{\left [ c_d^2 (t - t')^2 - h^2\right ] ^{\frac{1}{2}} } \mathrm{d} \boldsymbol{x}' \mathrm{d} t' \end{split} \end{equation}
where
$\Box ^2 = \nabla ^2 - \frac{1}{c_d} \frac{\partial ^2}{\partial t^2}$
is the wave operator,
$\nabla ^2 = \Sigma \frac{\partial ^2}{\partial x_i^2}$
is the Laplacian operator,
$c_d$
is the propagation speed of the sound wave in the material of the end actuator,
$\rho$
is the density of the material of the end actuator, and
$g(\boldsymbol{x})$
represents the external disturbance, which is the body force potential exerted on the end actuator at point
$\boldsymbol{x}$
.
$\phi$
represents the strain potential caused by the body force potential
$g$
, and the form of this general solution depends on the form of the external disturbance
$g$
.
3. Determination of tension distribution coefficients in flexible cable
The algorithm presented in this paper targets 3-DOF CDPRs utilizing four cables, for which an elastic dynamic mathematical model of the end effector is established below. For a planar three-degree-of-freedom CDPR, it usually has the structure shown in Figure 1.
Schematic diagram of the
$i$
th pulley-cable-end actuator structure.
Assuming that one of the forces acting on the end actuator is
$ \boldsymbol{F} = f(t) \delta (\boldsymbol{x} - \boldsymbol{p}) \boldsymbol{i}_1$
, where
$\delta (\boldsymbol{x})$
is a two-dimensional Dirac function. According to vector field decomposition,
$\boldsymbol{F}$
can be decomposed as follows:
\begin{align} & g = \nabla \cdot \boldsymbol{W}, \quad G = - \nabla \times \boldsymbol{W} \end{align}
\begin{align} & f = \nabla g + \nabla \times \boldsymbol{G} = \nabla ^2 \boldsymbol{W} \end{align}
where
$g$
represents the volumetric expansion part in the displacement potential and
$\boldsymbol{G}$
represents the rotational part in the displacement potential, which is not considered in this paper. Eq. (4) is the Poisson equation. According to the theory of partial differential equations, its fundamental solution can be expressed in the following form:
\begin{equation} \boldsymbol{W} (t) = \int _{V_x} \Phi (\boldsymbol{x} - \boldsymbol{x}^{\prime})\,f \, \mathrm{d} V \end{equation}
where:
\begin{equation} \varPhi (\boldsymbol{x}) = \left \{ \begin{split} &- \frac{1}{2 \pi } \ln{|\boldsymbol{x}|} &,n = 2\\ & \frac{\sqrt{\pi ^n}}{\alpha (n) \Gamma (\frac{n}{2}+ 1)} \frac{1}{|\boldsymbol{x}|} &,n\gt 2 \\ \end{split} \right . \end{equation}
For
$n=2$
, we have the following:
\begin{equation} \boldsymbol{W} (\boldsymbol{x},t) = -\frac{f(t)}{2 \pi } \ln{|\boldsymbol{x} - \boldsymbol{p}|} \boldsymbol{i}_1 \end{equation}
\begin{equation} g(\boldsymbol{x},t) = \nabla \cdot \boldsymbol{W} = -\frac{f(t)}{2 \pi } \frac{1}{r} \frac{\partial r}{\partial x_1} \end{equation}
\begin{equation} [g(\boldsymbol{x}^{\prime}, t)] = - \frac{f(t - \frac{\rho }{c_d})}{2 \pi \rho c_d ^2} \frac{\partial }{\partial x_1} (\!\ln{R}) \end{equation}
where
$r = \sqrt{(x_1 - p_1)^2 + (x_2 - p_2)^2}$
,
$R = \sqrt{(x^{\prime}_1 - p_1)^2 + (x^{\prime}_2 - p_2)^2}$
,
$\boldsymbol{i}_1$
represents the unit tensor or unit vector of
$\boldsymbol{x}_1$
, [] represents the translation operator, and
$[g(\boldsymbol{x}^{\prime}, t)]$
denotes the translation of the function
$g(\boldsymbol{x}^{\prime}, t)$
along the direction of
$\boldsymbol{x}_i$
. Substituting Eq. (2) yields:
\begin{equation} \begin{split} &\varPhi (\boldsymbol{x}, t;\; \boldsymbol{p}_i) = \frac{c_d}{2 \pi } \int _{V_d} \frac{ [g] }{ [c_d ^2 (t - t^{\prime}) ^2 - \rho ^2] ^{\frac{1}{2}} } \mathrm{d} V _{x^{\prime}} \\ &\quad \quad \quad = \frac{1}{4 \pi ^2 c_d \rho } \cdot \int \frac{f(t - t^{\prime} - \frac{h}{c_d})}{\sqrt{c_d ^2 (t - t^{\prime})^2 - R^2}} \frac{\partial }{\partial x^{\prime}_1} (\!\ln{R}) \mathrm{d} x^{\prime}_1 \mathrm{d} x^{\prime}_2 \mathrm{d} t^{\prime} \end{split} \end{equation}
On the circle determined by
$t'$
and
$h$
,
$\frac{f(t - t') - \frac{h}{c_d}}{\sqrt{c_d^2 (t - t')^2 - R^2}}$
is a constant value. Therefore, Eq. (10) can be rewritten as follows:
\begin{equation} \varPhi (\boldsymbol{x}, t;\, \boldsymbol{p}_i) = \frac{1}{4 \pi ^2 c_d \rho } \int \frac{f(t - t')}{\sqrt{c_d^2 (t - t')^2 - h^2}} \mathrm{d} h \cdot \int \frac{\partial }{\partial x^{\prime}_1} (\!\ln{R}) \mathrm{d} S \end{equation}
That is, first integrate along the circle determined by
$h$
and
$t'$
, and then integrate over
$h$
and
$t'$
.
Here,
$h = c_d (t - t')$
is a variable dependent only on time,
$\frac{\partial }{\partial x^{\prime}_1} (\!\ln{R})$
is an odd function with respect to
$x^{\prime}_1$
, and the integration domain
$S$
is a closed circular contour symmetric about
$x_1$
. Therefore, when
$\boldsymbol{p}_i$
is inside
$S$
, the second half of the integral is 0; otherwise, it is
$2 \pi R^2$
. Since the upper limit of the integration variable
$h$
can be replaced by
$r$
, and
$h = c_d (t - t')$
, Eq. (11) can be transformed to the following form:
\begin{equation} \varPhi (\boldsymbol{x}, t; \boldsymbol{p}_i) = \frac{1}{2 \pi c_d \rho } \frac{\partial }{\partial x_1} \ln{r} \int \frac{ f(t - t' - \frac{h}{c_d})}{\sqrt{c_d^2 (t - t')^2 - h^2}} \mathrm{d} t' \mathrm{d} h \end{equation}
Now let’s define the influence factor
$\varphi _i$
:
\begin{equation} \varphi _i (t) = \int _S \varPhi \mathrm{d} x_1 \mathrm{d} x_2 \end{equation}
Here,
$S$
represents the equivalent two-dimensional shape of the end effector. This influence factor represents the effect of applying force
$f$
at position
$p_i$
on the end effector.
Let’s define the influence coefficient
$\lambda _i$
:
\begin{equation} \lambda _i = \varphi _i \left ( \sum _i \varphi _i \right )^{-1} \end{equation}
Now, from the perspective of rigid body mechanics, let’s establish the dynamic equation of the end effector:
\begin{equation} \sum \boldsymbol{T}_i = m \sum \ddot{\boldsymbol{x}}_i = m \ddot{\boldsymbol{p}}_d \end{equation}
Here,
$T_i$
represents the tension of the
$i$
th cable.
For CDPRs, a commonly used approach to drive the end effector is by employing redundant cables, which exceed the degrees of freedom of the end effector. However, according to the theory of analytical mechanics, a planar rigid body can have at most three independent dynamic equations, while a planar CDPR with three degrees of freedom requires a minimum of four cables. In other words, there is no one-to-one mapping between pushing from left to right and pushing from right to left as shown in Eq. (15). To address this issue, numerous scholars have proposed different solutions. These methods typically involve iterative and computationally complex calculations. Moreover, when applied to practical control, they are limited by hardware performance and struggle to meet real-time requirements. In this paper, an attempt is made to define the influence of each force on the end effector from the perspective of elastic dynamics in order to calculate all driving forces. Consequently, the following formulas are provided:
\begin{equation} \begin{split} &\sum _i \boldsymbol{T} _i = m \sum _i \ddot{\boldsymbol{p}} _i = m \ddot{\boldsymbol{p}} _d \\ &\boldsymbol{p} _i = \lambda _i \boldsymbol{p} _d \end{split} \end{equation}
Alternatively, it can be written as:
\begin{equation} \boldsymbol{T} _i = \lambda _i m \boldsymbol{p} _d \end{equation}
Since
$\phi$
is calculated based on the components of externally applied forces, Eq. (17) needs to be expressed in component form as follows:
\begin{equation} T _{ij} = \lambda _{ij} mp _{dj} \end{equation}
Here,
$T_{ij}$
represents the
$j$
th component of the
$i$
th force
$T_i$
(distinct from the tensor). Considering that for CDPRs, the number of independent forces is equal to the number of independent cables, out of these eight components of
$T_{ij}$
, only four are independent. The relationship between these eight components will be established below.
It is not difficult to find that by substituting Eqs. (13) and (14) into Eq. (18), we obtain an integral equation about
$T_{ij}$
. To reduce the length, let
$\xi = \sqrt{c_d^2(t - t')^2 - h^2}$
as follows:
\begin{equation} T_{ij} = \frac{\int \left (\frac{\partial }{\partial x_1} \ln{r} \int \frac{[T_{ij}]}{\xi } \mathrm{d} t' \mathrm{d} h\right ) \mathrm{d} \Omega }{\sum _i \int \left (\frac{\partial }{\partial x_1} \ln{r} \int \frac{[T_{ij}]}{\xi } \mathrm{d} t' \mathrm{d} h\right ) \mathrm{d} \Omega } m \ddot{\boldsymbol{p}}_{dj} \end{equation}
Since Eq. (19) is complex and difficult to solve, it is simplified as follows:
\begin{equation} T_{i1} = \varphi _{i1} m \ddot{p}_{\mathrm{d}1} \end{equation}
Making a variable substitution for
$\varphi _{ij}$
:
$u = c_d(t - t')$
, where
$u \in (c_d t, 0)$
, we obtain:
\begin{equation} \varphi _{ij} = \int \frac{\partial }{\partial x_1} \ln{r_i} \int \frac{[T_{ij}]}{\sqrt{u^2 - h^2}} \mathrm{d} h \mathrm{d} u \mathrm{d} \Omega \end{equation}
Continuing the variable substitution by letting
$h = u \sin{\gamma }$
and
$dh = u \cos{\gamma }\mathrm{d} \gamma$
, where
$\gamma \in \left ({-}\frac{\pi }{2}, \frac{\pi }{2}\right )$
, we have:
\begin{equation} \varphi _{ij} = \int \frac{\partial }{\partial x_1} \ln{r_i} \int [T_{ij}] \mathrm{d} \gamma \mathrm{d} u \mathrm{d} \Omega \end{equation}
At this point,
$[T_{ij}] = T_{ij}\left (\frac{u - u\sin{\gamma }}{c_d}\right )$
. It can be verified that
$\int [T_{ij}] \mathrm{d} \gamma \mathrm{d} u$
is a function only dependent on time
$t$
, denoted as
$G(t)$
. Therefore, we obtain:
\begin{equation} \varphi _{ij} = kG(t) \end{equation}
Clearly,
$k = \int \frac{\partial }{\partial x_1} \ln{r_i} \mathrm{d} \Omega$
is a physical parameter related to the shape of the end effector, which is constant for a given end effector. Substituting Eq. (23) back into Eq. (22), we have:
\begin{equation} kG_{ij}(t) = k \int _{c_dt}^0 \mathrm{d} u \int _{-\frac{\pi }{2}}^{\frac{\pi }{2}} k[G_{ij}] m [\ddot{p}_{dj}] \mathrm{d} \gamma \end{equation}
Taking the derivative of Eq. (24) with respect to
$t$
on both sides, we obtain:
\begin{equation} \dot{G}_{ij}(t) = -C_{\mathrm{P}} \int _{-\frac{\pi }{2}}^{\frac{\pi }{2}} [G_{ij}][\ddot{p}_{dj}] \mathrm{d} \gamma \end{equation}
In this case,
$C_{\mathrm{P}} = m\,k\,c_d$
, and
$[G_{ij}] = G(t(1 - \sin{\gamma }))$
. Continuing with the variable substitution in Eq. (25), let
$v = t(1 - \sin{\gamma })$
. Then
$\gamma = \arcsin{(1 - \frac{v}{t})}$
, and
$\mathrm{d} \gamma = \frac{1}{\sqrt{1 - \left ( 1 - \frac{v}{t} \right ) ^2}} \mathrm{d} v$
yields:
\begin{equation} \dot{G}_{ij} (t) = -C_{\mathrm{P}} \int _{2t} ^0 \frac{G_{ij} (v) \ddot{p} _{dj} (v)}{\sqrt{1 - \left ( 1 - \frac{v}{t} \right ) ^2}} \mathrm{d} v \end{equation}
Eq. (26) is complex and difficult to solve. However, for engineering problems, it is often sufficient to find a particular solution that satisfies the engineering conditions. Due to the appealing properties of polynomials, this paper attempts to find a set of solutions that satisfy the aforementioned equation from the polynomial space. This paper provides the following set of solutions: when
$G_{ij} (t)$
can be decomposed into the product of two rational polynomials,
$g_1(t)$
and
$g_2(t)$
, where
$g_2(t)$
is a quadratic polynomial, and let
$\ddot{p}_{dj} = g_2(t) ^{-1}$
, then there must exist
$G_{ij}(t)$
and
$p_{dj} (t)$
that satisfy Eq. (26).
Proof: By performing a variable substitution in Eq. (26), let
$\sin{w} = 1 - \frac{v}{t}$
, then
$\mathrm{d} v = -t\,\cos{w}\,\mathrm{d} w$
, and
$v = t\,(1 - \sin{w})$
.
\begin{equation} \begin{split} &\dot{G} _{ij} (t) = C_{\mathrm{P}}\,t \int _{-\frac{\pi }{2}} ^{\frac{\pi }{2}} G_{ij}(t - t\sin{w}) \ddot{p} _{dj} (t - \sin{w}) \mathrm{d} w \\ &\frac{\mathrm{d}}{\mathrm{d} t} (g_1(t) g_2(t)) = C_{\mathrm{P}}\,t \int _{-\frac{\pi }{2}} ^{\frac{\pi }{2}} g_1(t -t\sin{w}) \mathrm{d} w \end{split} \end{equation}
Here,
$g_1(t - t\,\sin{w})$
is a polynomial in terms of
$t - t\,\sin{w}$
, and assume the
$i$
th term of
$g_1$
is
$a_i (t - t\,\sin{w})^i = a_i t^i \mathrm{C}_i ^m({-}\sin{w})^{i-m}$
. Therefore:
\begin{align} &a_i t^i \int _{-\frac{\pi }{2}} ^{\frac{\pi }{2}} \mathrm{C}_i ^m({-}\sin{w})^{i-m} \mathrm{d} w = \nonumber\\ & \begin{cases} 0,\quad \text{when $i$ and $m$ have different parity} \\ \frac{i - m -1}{i - m} \ldots \frac{1}{2} \cdot \pi \, \mathrm{C}_i ^m a_i t^i,\text{otherwise} \end{cases} \end{align}
It can be seen that after integrating the right side of Eq. (27), it remains a polynomial in terms of
$t$
, with a degree one higher than the original
$g_1(t)$
. Therefore, there should exist a polynomial
$G_{ij}$
of degree two higher than
$g_1$
that satisfies Eq. (27). In other words, let
$g_2$
be a quadratic polynomial in
$t$
, and Eq. (27) can be satisfied. Since the proof for higher-degree polynomials is not the focus of this paper, we only prove the existence of
$g_1$
in quadratic polynomial form.
Let
$g_1 = b_2 t^2 + b_1 t + b_0$
and
$g_2 = (\ddot{p}_{d1} (t))^{-1} = c_1 t^2 + c_2 t + c_3$
. Based on the equality of corresponding coefficients, the following system of equations can be obtained:
\begin{equation} \mathbf{A}\,b = 0 \end{equation}
where
\begin{align} \mathbf{A}=&\left [ \begin{array}{c@{\quad}c@{\quad}c} 0& 4c_3& 4c_2-\text{C}_{\mathrm{P}}\\ 3c_3-C_{\text{p}}& 3c_2& 3c_1\\ c_2& c_1& c_0\\ c_1& c_0& 0 \end{array} \right ] \nonumber\\ \boldsymbol{b}=&\left [ \begin{array}{c@{\quad}c@{\quad}c} b_0& b_1& b_2 \end{array} \right ] \end{align}
It can be seen that in order for the above equations to have a solution, the rank of the matrix must be less than 3. When the relationships between the elements of
$\boldsymbol{c}=\left [ \begin{matrix} c_0& c_1& c_2\\ \end{matrix} \right ]$
satisfy the conditions in the appendix, the rank of matrix
$\mathbf{A}$
is 2. That is, when the values of
$b_1$
,
$b_2$
,
$b_3$
,
$c_1$
,
$c_2$
, and
$c_3$
satisfy the conditions in the appendix, the assumption holds and the proof is complete.
Through the above proof, it can be found that the
$G_{ij}(t)$
determined by
$p_{d1}$
is not unique. This is because the matrix
$\mathbf{A}$
does not have full rank, and any linear combination of its null vectors satisfies the equation set (29). From Eq. (27), it can be easily seen that any multiple
$\beta$
of the null vector can be absorbed into
$C_{\mathrm{P}}$
. Therefore, the exact value of
$C_{\mathrm{P}}$
is not important in the derivation process of the preceding formulas, as there always exists a suitable
$\beta$
to make
$C_{\mathrm{P}}$
take any value.
3.1. Approximate calculation method for cable tension direction
Illustration of tangent point and approximation point.
After determining the tension components in either the horizontal or vertical orientation of the cable, the corresponding tension vector should be derived from the cable’s direction vector. An essential element in this process is the Jacobian matrix. Owing to the complexity of its rigorous calculation, this paper presents a simplified method for computing the Jacobian matrix.
As shown in Figure 2, under normal circumstances, the actual contact point slides within the dashed region of the
$\frac{1}{4}$
arc. Calculating the Jacobian matrix directly using the real position of the contact point
$\boldsymbol{p}_{\mathrm{r}}$
can be cumbersome. However, if we allow for some error, we can simplify the calculation of the cable tension direction and the subsequent Jacobian matrix by artificially fixing the contact point at position
$\boldsymbol{p}_{45}$
on the arc, which forms a
$45^\circ$
angle with the
$x$
-axis. Since the cable direction vectors in the Jacobian matrix are normalized, this paper evaluates the algorithm’s error using angles. Let
$\boldsymbol{p}_{\mathrm{r}}$
be the actual contact point position,
$\boldsymbol{p}_{45}$
be the artificially defined contact point position, and
$\boldsymbol{p}_{\mathrm{e}}$
be the end effector position. Based on geometric relationships, the maximum error in
$\theta$
occurs at the position shown in the right figure of Figure 2. Assuming counterclockwise deviation is positive and clockwise deviation is negative, it can be inferred that the deviation lies within the range shown in Eq. (31).
\begin{equation} \begin{split} & \theta _e ^{-} = 0 \\ & \theta _e ^{+} = \arcsin{\left(\frac{R \left (1 - \frac{\sqrt{2}}{2}\right )}{\|\boldsymbol{p}_{\mathrm{e},2} - \boldsymbol{p}_{45}\|}\right)} \end{split} \end{equation}
According to Eq. (31), the maximum deviation value is independent of the pulley’s position. If we place the pulley at
$({-}R, -R)$
, then
$\boldsymbol{p}_{45}$
will be at the origin. In this case,
$\theta ^+ _e = R(1-\frac{\sqrt{2}}{2})/|\boldsymbol{p}_{e2}|$
, resulting in Figure 3. When the ratio is 10.2, the error is 1.6 degrees, and when the ratio is 15, the error is 1.1 degrees.
Variation of
$\theta ^+ _e$
with
$R/\|\boldsymbol{p}_{\mathrm{e},2}\|$
.
Let
$\boldsymbol{u}_i = u_{ij}\,(i = x,y;\,j=1,2)$
be the direction vectors of the cable tension. We have:
\begin{equation} \boldsymbol{u}_i = \frac{\boldsymbol{p}_{\mathrm{e}} - \boldsymbol{p}_{45}}{\|\boldsymbol{p}_{\mathrm{e}} - \boldsymbol{p}_{45}\|} \end{equation}
Therefore, the relationship between the magnitude and direction of the cable tension can be written as follows:
\begin{equation} T_{ij} = \| \boldsymbol{T}\| _{i} u_{ij}\,(i=x,y;\ j=1,2) \end{equation}
3.2. The algorithm process
The specific process of the cable tension distribution algorithm introduced in this paper is depicted in Figure 4. It necessitates an understanding of the global desired trajectory and its second derivative to compute the coefficients for the time-varying tension distribution of the cables. This calculation can be considered pre-processing of the input signal, not part of the control strategy loop.
Unlike the cable tension distribution algorithms referenced in the introduction, which fully calculate the Jacobian matrix within the control strategy loop, this algorithm’s optimization branch employs an iterative solution based on optimization theory, while the convex geometry branch uses the geometric attributes of shapes such as vertices, edges, and centroids for its calculations. These processes occur within the control strategy loop, endowing the proposed algorithm with a natural real-time performance advantage.
Flowchart of the tension distribution algorithm.
4. Numerical Example
In this section, the aforementioned formulas are used to calculate a practical problem, and the influence of
$\beta$
on cable tension is discussed. The values of the main parameters in this paper are shown in Table I.
Values of the main parameters
4.1. Influence of
$c_2$
and
$C_{\mathrm{P}}$
on end trajectory
By setting the rank of matrix
$\mathbf{A}$
to 2, it can be determined that the parameters of
$\boldsymbol{c}$
can be expressed as functions of
$c_2$
and
$C_{\mathrm{P}}$
. Therefore, their values directly affect the value of the end trajectory. Finally, it is determined that
$c_2 = 2$
and
$C_{\mathrm{P}} = 4$
.
Influence of
$C_{\mathrm{P}}$
on the trajectory.
Influence of
$c_2$
on the trajectory.
4.2. Influence of
$\beta$
values on cable tension
Based on the values of the coefficients mentioned above, the specific values of each quantity can be calculated as follows (Figs. 5 and 6):
\begin{align} \mathbf{A} = \left [ \begin{array}{c@{\quad}c@{\quad}c} 0.000 & 0.533 & 12.858 \\[3pt] -1.600 & 12.000 & -82.654 \\[3pt] 4.000 & -27.551 & 265.670 \\[3pt] -27.551& 265.670 & 0 \end{array} \right ] \end{align}
\begin{align} \boldsymbol{b} = \boldsymbol{null}(\mathbf{A}) = \left [ \begin{array}{c@{\quad}c@{\quad}c} 0.995 & 0.103 & -0.004 \end{array} \right ] \end{align}
So that:
\begin{equation} \begin{split} g_{12}(t) =&\ g_{22}(t) \\ =&\ \ddot{p}^{-1} _{\mathrm{d}1} = 0.27\,t^3 +2.00\,t^2 +4.88\,t+4.74 \end{split} \end{equation}
\begin{equation} \begin{split} g_{11}(t) &= \beta _1 ({-}0.004\,t^2 + 0.103\,t+0.995) \\ g_{22}(t) &= \beta _2 ({-}0.004\,t^2 - 0.103\,t+0.995) \end{split} \end{equation}
Therefore, by substituting Eqs. (36), (37) into (20), (23), and (27), we can obtain the equation
$\lambda _{i1} = \lambda _{i1}(\beta _1, \beta _2, t)$
. This is a surface plot in four-dimensional space, making it difficult to obtain intuitive information. Hence, this paper uses
$\lambda _{i1,\mathrm{max}}$
to reflect the influence of
$\beta$
on
$\lambda _{ij}$
. Here,
$\lambda _{i1,\mathrm{max}}$
represents the extremum in the time direction. For the specific derivation process, please refer to the Appendix C. The three-dimensional surface plot can be referred to in Figure 7. It is easy to see that different rope tensions will directly affect the maximum value of the rope tension in the trajectory planning. Generally, the greater the rope tension, the stronger the anti-interference ability of the CDPR, so the rope tension should be as large as possible within the allowable output torque range of the motor. Considering all factors, the values of
$\beta _1$
and
$\beta _2$
can be chosen as shown in Table I.
Influence of
$\beta$
on
$\lambda _{11}$
.
4.3. Calculation of
$\theta _d$
When the cable tensions and end displacement are determined, the angle of the end effector is no longer arbitrary. This is because the angular displacement of the end effector of a CDPR (cable-driven parallel robot) is coupled with its displacement, requiring the calculation of the end effector’s angular displacement based on the cable tensions and end trajectory. By applying the moment balance equation to the end effector, we have:
\begin{equation} \begin{split} m\,\ddot{\theta }_d &= \sum _i f_i(\boldsymbol{r}_i \times \boldsymbol{u}_i) \\ &= \sum _i f_i(r_{i1}\,u_{i2} - r_{i2}\,u_{i1}) \end{split} \end{equation}
Here,
$\boldsymbol{r}_i$
and
$\boldsymbol{u}_i$
represent the position vector and unit vector of cable
$i$
, respectively. Thus, it can be known that:
\begin{equation} \theta _d = \int _0 ^t \left ( \int _0 ^t \ddot{\theta }_d \mathrm{d}t \right ) \mathrm{d} t \end{equation}
For an arbitrarily complex trajectory, its
$\theta _d$
generally cannot be analytically solved. The reason is that
$\boldsymbol{u}_i$
in Eq. (38) is a function of
$\theta$
. Therefore, in the subsequent simulation, numerical integration will be used for Eq. (39) (Figure 8).
4.4. Simulation analysis
Based on the previous selection, the expression for
$p_{d1}$
can be reverse calculated as follows:
\begin{equation} \begin{split} \boldsymbol{p} _{\mathrm{d}1} &= \left ( p_{\mathrm{d}1},\,p_{\mathrm{d}2},\,p_{\mathrm{d}3} \right )\\ p_{\mathrm{d}1} &= 0.53\,{\left (t+4.08\right )}\,\log \left (t+4.08\right ) \\ &-0.62\,\log \left ({{\left (t+1.71\right )}}^2 +1.428\right ) \\ & -0.26\,{\left (t+1.71\right )}\,\log \left ({{\left (t+1.71\right )}}^2 +1.43\right ) \\ & -0.63\,{\left (t+1.71\right )}\,\mathrm{atan}\left (0.84\,t+1.43\right ) \\ & +1.04\,\mathrm{atan}\left (0.84\,t+1.43\right ) - 1.2\,t - 3.5 \\ p_{\mathrm{d}2} &= p_{\mathrm{d}3} = 0 \end{split} \end{equation}
The values of
$\lambda _{i1}$
at this time are as follows:
\begin{equation} \begin{split} &\lambda _{11} = \lambda _{41} = \frac{-2.16\,t^2 -3.47\,t+3.37}{3.14\,t^2 -12.00\,t+11.64} \\ &\lambda _{21} = \lambda _{31} = \frac{4.77\,t^2 -7.68\,t+7.45}{2.82\,t^2 -10.80\,t+10.48} \end{split} \end{equation}
Substituting into Eq. (20), we can obtain the values of
$T_{11}$
and
$T_{41}$
:
\begin{equation} \begin{split} &T_{11} = T_{41} = \frac{345.37\,t^2 -555.87\,t+539.30}{T_{\mathrm{den}}} \\ &T_{21} = T_{31} = \frac{847.13\,t^2-1363.46\,t+ 1322.78}{T_{\mathrm{den}}} \\ &T_{\mathrm{den}} = 8.47\,t^5+30.33\,t^4- 55.46\,t^3 \\ &\qquad \quad -203.99\,t^2-6.72\,t+ 551.64 \end{split} \end{equation}
Substituting Eq. (42) into Eq. (33), the tension
$|\boldsymbol{T}_i|$
of the cable can be obtained. Furthermore, dynamic simulations yield the results shown in Figure 9(a). Figure 9(b) displays the simulated trajectory and error obtained using the Newton Slack Any P method [Reference Ueland, Sauder and Skjetne30]. Here,
$p_{\mathrm{ci}}$
represents the simulated trajectory in the
$i$
th direction, and
$p_{\mathrm{eri}}$
represents the error between the simulated trajectory and the desired trajectory
$(i=x,y)$
.
Tension of the cable in the
$x$
direction.
(a) Simulated trajectory and error in the
$i$
th direction using the proposed method. (b) Simulated trajectory and error in the
$i$
th direction using the Newton Slack Any P method
$(i=1,2)$
.
From Figure 8, it can be observed that the cable tension calculated using the proposed method is continuous. Figure 9 demonstrates that the results obtained from dynamic simulations using the proposed tension distribution method align well with the desired trajectory. Running the simulation trajectory 20 times and taking the average runtime yields an average time of 31.44 s for the simulation model using the proposed cable tension distribution method, while the average runtime for the simulation program using the Newton Slack Any P method is 35.37 s.
5. Conclusion
From the wave equation, this paper defines the influence factor
$\phi _i$
of the
$i$
th cable on the end effector and defines the influence coefficient
$\lambda _{ij}$
to allocate the cable tension. A mathematical model is established for
$\lambda _{ij}$
and
$\ddot{p}_{d1}$
, and it is proven that it has a solution in at least cubic polynomial form. This solution is a rational polynomial fraction form, which is computationally simple and can better meet real-time force allocation requirements. Finally, a specific numerical example is provided to verify the correctness of the algorithm and demonstrate its real-time capability. A simple comparison is made with the Newton Slack Any P method proposed by scholar Einar Ueland. However, when using this algorithm, if the position of the end effector is specified, its angle cannot be freely determined. Part of the reason is that the solution provided by this algorithm is not the solution in the entire space, but a curve in the solution space. Another reason is that the position and angle of CDPRs are strongly coupled, and specifying the displacement and angular displacement of the end effector does not necessarily result in a force helix curve within the feasible space of CDPRs.
Our team’s future work includes the following: (1) Planning to change
$\lambda _{ij}$
to expand the range of trajectories applicable to this algorithm, aiming to make this algorithm usable for any curve within the feasible space of force helixes in CDPRs. (2) Planning to introduce a friction model to more comprehensively describe and predict the behavior of cables in actual operations. (3) Expanding the range of applications for the algorithm, not limited to four-cable 3-DOF CDPRs, extending to spatial CDPRs.
Author contributions
Da Song, Ming Lu, and Lei Zhao completed the main part of the theoretical derivation. Zhichao Sun and Haochen Wang participated in the theoretical simulation and translation of the paper. Lixun Zhang provides special guidance in the theoretical part of this paper.
Financial support
Supported by the National Natural Science Foundation of China (Grant No. 61773007) and the Jilin Provincial Natural Science Foundation of China (Grant No. YDZJ202301ZYTS272)
Competing interests
The authors declare no conflicts of interest exist.
Ethical approval
Not applicable.
Appendix A. Condition for Rank-2 of
$\mathbf{A}$
There are multiple methods to determine the rank-2 condition of
$\mathbf{A}$
, and this paper provides the following method:
\begin{equation} \begin{split} a_{c1} &=-\frac{16\,{C_{\mathrm{P}}}^2 }{75} \\ a_{c2} &= -\frac{4\,C\mathrm{p}\,{{\left (4\,c_2 -\frac{\pi \,C_{\mathrm{P}}}{2}\right )}}^2 }{5} \\ &\quad -\frac{8\,C\mathrm{p}\,{\left (4\,{c_2 }^2 -\frac{\pi \,C_{\mathrm{P}}\,c_2 }{2}\right )}}{5} \\ &\quad +\frac{8\,C_{\mathrm{P}}\,c_2 \,{\left (4\,c_2 -\frac{\pi \,C_{\mathrm{P}}}{2}\right )}}{5} \\ a_{c3} &= 6\,c_2 \,{\left (4\,c_2 -\frac{\pi \,C_{\mathrm{P}}}{2}\right )}\,{\left (4\,{c_2 }^2 -\frac{\pi \,C_{\mathrm{P}}\,c_2 }{2}\right )} \\ &\quad -3\,{{\left (4\,{c_2 }^2 -\frac{\pi \,C_{\mathrm{P}}\,c_2 }{2}\right )}}^2 \\ c_1 &= \frac{-a_{c2} +- \sqrt{a_{c2} ^2 - 4\,a_{c1}\,a_{c3}}}{2\,a_{c1}} \\ c_0 &= \frac{4\,c_1 ^2\,c_2 - \frac{\pi }{2}c_1 ^2\,C_{\mathrm{P}}}{4\,c_1\,c_3 + 4\,c_2 ^2 - \frac{\pi }{2}C_{\mathrm{P}}} \end{split} \end{equation}
Where
$c_2$
and
$C_{\mathrm{P}}$
are arbitrary real numbers, and
$c_0$
and
$c_1$
determined by Eq. (43) make the rank of
$\mathbf{A}$
equal to 2.
Appendix B. Integral form of
$\ddot{p}_{\mathrm{d}1}$
in the real number domain
For the integration of
$\ddot{p}_{\mathrm{d}1}$
, the most complex case arises when
$\ddot{p}_{\mathrm{d}1}^{-1}$
has only one intersection with the
$x$
-axis (which is the case provided in this paper). In this situation, it has a pair of complex conjugate roots, and the integration of
$\dot{p}_{\mathrm{d}1}$
obtained by conventional methods is a complex solution. However, from a graphical perspective, the area enclosed by
$\ddot{p}_{\mathrm{d}1}$
,
$t=0$
, and the
$t$
-axis is a real number (a real function). From a mathematical point of view, the solution in the complex number domain (complex function solution) can always be simplified into a complex number (complex function) with a nonzero real part and a zero imaginary part. However, this is often difficult or even impossible to achieve. Therefore, in this paper, an integration method is proposed to avoid using complex numbers.
Assume that when
$\ddot{p}_{\mathrm{d}1}(t)$
has one real solution
$t_0$
and a pair of complex conjugate solutions
$t_r + t_i i$
and
$t_r - t_i i$
, the following equation holds:
\begin{equation} \ddot{p}_{\mathrm{d}1}(t) = \frac{1}{c_3 (t - t_0)((t-t_r)^2 + t_i^2)} \end{equation}
Assume:
\begin{equation} \ddot{p}_{\mathrm{d}1}(t) = \frac{1}{c_3}\left (\frac{e_1}{t - t_0} + \frac{e_2 t + e_3}{(t-t_r)^2 + t_i^2}\right ) \end{equation}
By equating the coefficients in Eqs. (44) and (45), the following system of Eq. (46) is obtained:
\begin{equation} \mathbf{A}_{\mathrm{p}} \,\boldsymbol{e} = [\begin{array}{c@{\quad}c@{\quad}c} 0 & 0 & 1 \end{array}]^{\mathrm{T}} \end{equation}
Where:
\begin{align} & \mathbf{A}_{\mathrm{p}} = \left [\begin{array}{c@{\quad}c@{\quad}c} 1 & 1 & 0\\ -2\,t_r & -t_0 & 1\\{t_i }^2 +{t_r }^2 & 0 & -t_0 \end{array}\right ] \end{align}
\begin{align} & \boldsymbol{e} = \left [\begin{array}{c@{\quad}c@{\quad}c} e_1 & e_2 & e_3 \end{array}\right ]^{\mathrm{T}} \end{align}
Solving Eq. (46) yields:
\begin{equation} \boldsymbol{e} = \mathbf{A}_p^{-1}\left [\begin{array}{c} 0 \\ 0 \\ 1 \end{array}\right ] = \left [\begin{array}{c} \frac{1}{{t_i }^2 +{t_r }^2 +2\,t_0 \,t_r -t_0 }\\[3pt] -\frac{1}{{t_i }^2 +{t_r }^2 +2\,t_0 \,t_r -t_0 }\\[3pt] \frac{2\,t_r}{{t_i }^2 +{t_r }^2 +2\,t_0 \,t_r -t_0 } \end{array}\right ] \end{equation}
At this point, the difficult-to-integrate part in Eq. (45) is reduced to the second half, given by:
\begin{equation} \frac{e_2 t + e_3}{(t-t_r)^2 + t_i ^2} = \frac{e_2}{2}\left (\frac{2\,t-2\,t_r}{(t-t_r)^2 + t_i ^2} + \frac{2\,e_3/e_2 + 2\,t_r}{(t-t_r)^2 + t_i ^2}\right ) \end{equation}
At this point, it can be easily obtained that:
\begin{equation} \begin{split} \dot{p}_{\mathrm{d}1} =& \int _0 ^t \ddot{p}_{\mathrm{d}1} (\tau ) \mathrm{d} \tau \\ =& \frac{1}{c_3} \int _0 ^t \frac{e_1}{\tau - t_0} \mathrm{d}\tau + \frac{e_2}{2\,c_3} \int _0 ^t \frac{2\,\tau - 2\,t_r}{(\tau -t_r)^2 + t_i ^2} \mathrm{d} \tau \\ &+\frac{e_3+e_2\,t_r}{c_3} \int _0 ^t \frac{1}{(\tau -t_r)^2 + t_i ^2} \mathrm{d} \tau \end{split} \end{equation}
After rearranging, we obtain:
\begin{equation} \begin{split} \dot{p}_{\mathrm{d}1} =& \frac{e_1}{c_3} \ln{(t - t_0)} + \frac{e_2}{2\,c_3} \ln{((t-t_r)^2 + t_i ^2)} \\ &+ \frac{e_3+e_2\,t_r}{c_3\,t_i} \mathrm{atan}{\frac{t-t_r}{t_i}} \end{split} \end{equation}
Continuing the integration, we have:
\begin{equation} \begin{split} p_{\mathrm{d}1}=&\frac{e_1}{c_3}\left ( -t+\left ( t-t_0 \right ) \ln \left ( t-t_0 \right ) \right )\\ &+\frac{e_2}{2\,c_3}\left ( -2\,t+2\,t_i\,\text{atan}\left ( \frac{t-t_r}{t_i} \right ) \right .\\ &\left . +\log \left ( \left ( t-t_r \right ) ^2+t_i^2 \right ) \,\left ( t-t_r \right ) \right )\\ &+\frac{e_3+e_2\,t_r}{c_3\,t_i}\left ( \left ( t-t_r \right ) \text{atan}\frac{\left ( t-t_r \right )}{t_i} \right .\\ &\left . -\frac{t_r}{2}\ln \left ( \left ( t-t_r \right ) ^2+t_{i}^{2} \right ) \right )\\ \end{split} \end{equation}
Combining Eqs. (49) and (53), we obtain the real function integral solution of the integrand with a third-degree polynomial in the denominator, which has complex conjugate roots in the numerator.
Appendix C. Calculation method for
$\lambda _{i1,\mathrm{max}}$
Substituting
$g_{11}$
and
$g_{21}$
into Eqs. (23) and (14), we can obtain:
\begin{equation} \lambda _{i1} = \lambda _{i1} (b_0, b_1, b_2, \beta _1, \beta _2, t) \end{equation}
Taking the partial derivative of the above equation with respect to
$t$
, we have:
\begin{equation} \begin{split} \frac{\partial \lambda _{i}}{\partial t} &= \frac{\left (-2\,b_1 \,b_2 \,\beta _1 \,\beta _2 \right )\,t^2 +\left (-4\,b_0 \,b_2 \,\beta _1 \,\beta _2 \right )\,t}{d_4\,t^4 + d_3\,t^3 + d_2\,t^2 +d_1\,t + d_0 } \\ d_4 &= \left ({b_2 }^2 \,{\beta _1 }^2 +2\,{b_2 }^2 \,\beta _1 \,\beta _2 +{b_2 }^2 \,{\beta _2 }^2 \right ) \\ d_3 &= \left (2\,b_1 \,b_2 \,{\beta _1 }^2 -2\,b_1 \,b_2 \,{\beta _2 }^2 \right ) \\ d_2 &= \left ({b_1 }^2 \,{\beta _1 }^2 -2\,{b_1 }^2 \,\beta _1 \,\beta _2 +{b_1 }^2 \,{\beta _2 }^2 \right . \\ &\quad \left . +2\,b_0 \,b_2 \,{\beta _1 }^2 -2\,b_0 \,b_2 \,{\beta _2 }^2 \right ) \\ d_1 &= \left (2\,b_0 \,b_1 \,{\beta _1 }^2 -4\,b_0 \,b_1 \,\beta _1 \,\beta _2 +2\,b_0 \,b_1 \,{\beta _2 }^2 \right ) \\ d_0 &={b_0 }^2 \,{\beta _1 }^2 -2\,{b_0 }^2 \,\beta _1 \,\beta _2 +{b_0 }^2 \,{\beta _2 }^2 \end{split} \end{equation}
It is not difficult to observe that Eq. (55) has monotonic intervals with respect to
$t$
as
$\big(0, -\frac{b_0}{b_1}\big)$
and
$\big({-}\frac{b_0}{b_1}, +\infty \big)$
. Therefore, substituting
$t = -\frac{b_0}{b_1}$
along with
$b_0$
,
$b_1$
, and
$b_2$
back into Eq. (54), we obtain the three-dimensional surface plot of
$\lambda _{i1,\mathrm{max}} = \lambda _{i1} (\beta _1, \beta _2)$
, showing the maximum value on
$t$
varying with
$\beta _1$
and
$\beta _2$
.
