**4. NUMERICAL SOLUTION PROCEDURES AND VIRTUAL PROTOTYPE**

**4.1.1.1 Initial condition search**

an additional constraint equation for the length of a virtual link π·π·_{1}πΉπΉ_{1} form a set of
nonlinear, non-differential, algebraic equations. These equations are used to find the initial
conditions of the trajectories and to perform the required kinematics analysis. The nature
of these equations makes the direct search for initial conditions a challenging task and
requires a special numerical treatment. The differentiation of these equations with respect
to time yields a set of DAE that can be solved using advanced numerical methods (Newton-
Raphson) since it provides fast convergence to the solution.

**β’ Full set of constraint equations. A full representation of the constraint algebraic **
equations of the trajectories (ππ_{4}, ππ_{5}, ππ_{6}) of the front-end assembly, as well as the linear
displacements of the hoist and drag ropes (ππ_{7}, ππ_{8}) is given in a set of equations (4.1). The
symbolic abbreviations, βcβ and βsβ denote the trigonometric functions, sin and cos,
respectively. The linear displacements functions of the hoist and drag ropes are input into
the kinematics model and they take the form of the last two rows of equations (4.1).

76.60 + 1.715ππ4β 10.5ππ5β ππ6ππ7+ ππ8π π 4 = 0

56.89β + ππ4ππ8β 1.715π π 4β 10.5π π 5+ ππ7π π 6 = 0

β11986.81 + 220 ππ οΏ½37ππ_{180}+ ππ6οΏ½ ππ7β ππ72+ ππ82β 21 ππ[1.57β + ππ4+ ππ5]οΏ½2.94β + ππ82=0

β75 + 1.32π‘π‘ + ππ_{7} = 0

β75 β 2.54π‘π‘ + ππ8 = 0

To solve the kinematics and dynamics models of the dragline, all of the constraint
equations must be satisfied during the numerical integration. Any violation will force the
numerical solver to stop the integration, otherwise the numerical results maybe erroneous.
These constraint equations define the limits introduced by the machine geometry during
the digging phase and are then modified to meet the loaded bucket swinging motion. The
initial conditions of (ππ_{4}, ππ_{5}, ππ_{6}) can be resolved by switching these algebraic constraint
equations into a set of DAEs that accept initial values of the kinematics entities
(ππ_{4}, ππ_{5}, ππ_{6}, ππ_{7}, ππ_{8}). Equation (4.2) represents a new set of DAEs from the direct
differentiation of equation (4.1) with respect to the independent variable (t).

(c_{4}q_{8} β 1.715s_{4})qΜ_{4}+ 10.5s_{5}qΜ_{5}+ q_{7}s_{6}qΜ_{6}β c_{6}qΜ_{7}+ s_{4}qΜ_{8} = 0
β(1.715c_{4} + q_{8}s_{4})qΜ_{4}β 10.5c_{5}qΜ_{5}+ c_{6}q_{7}qΜ_{6}+ s_{6}qΜ_{7}+ c_{4}qΜ_{8} = 0
21 s[1.57 + q_{4}+ q_{5}]οΏ½2.94 + q2_{8} (qΜ_{4}+ qΜ_{5}) β 220 s οΏ½37Ο
180+ q6οΏ½ q7qΜ6
+220 c οΏ½37Ο
180+ q6οΏ½ qΜ7β 2q7qΜ7+ 2q8qΜ8β
21 c[1.57+q4+q5]q8qΜ8
οΏ½2.94+q_{8}2 = 0
1.32β + qΜ_{7} = 0
β2.54 + qΜ_{8} = 0
(4.2)
(4.1)

Equation (4.2) takes the form of equation (4.3), which is very common in constrained multibody mechanical systems.

Fππ(qi, qΜi, t) = 0 (i, j) = (1, . . ,5) (4.3)

Fππ represents the jth DAE, qi and qΜi represent the dependent variables or state variables

and the time derivative, respectively. An important step to follow when solving the DAEs
set in equation (4.2) is to select a good intital value of each dependent variable. This feature
defines the difference between the numerical solutions of the DAEs and ODEs. The latter
does not necessarily have a good initial value because no hidden constraints must be
satisfied at each integration step. A major difficulty in arriving at a quick numerical
solution has been attributed to the nature of the dependent variables. In other words, the
generalized coordinates do not have the same units and other entities (ππ_{7} and ππ_{8}) vary
signficantly with time.

Several numerical experiments have been carried out in Matlab and Mathematica to understand the role of every dependent variable on the numerical stability and convergence of the solution algorithm. If no initial values are given to the numerical alogrithm, NDSolve searches for initial values that yields zero DAEs residual. The numerical alogrithm, based on the Newton-Raphson method, achieved a good convergence at 1.52 seconds, but failed to advance the integration due to the presence of singularity at this time step. The trajectories of all ropes based on the initial values are shown in Figure 4.2 (a), with quick changes in their values at the onset of digging. These results show further enhancement of the numerical experiments and can be used as a starting point to construct new values.

Figure 4.2 (b) shows the variations in the lengths of the hoist and drag ropes, with
a linear behavior as the bucket penetrates the ground. It can be seen that the initial values
of the kinematics parameters, in radians, are ππ_{4,1}[0] = 10 β Pi 180β , ππ_{5,1}[0] = 30 β
Pi 180β , and ππ6,1[0] = β30 β Pi 180β , respectively, and the initial values of the linear

displacements of hoist and drag ropes are chosen as 75 m. These linear displacedments are found to reduce the numerical intgeration error. From an operational viewpoint, the rapid change of the angular displacement of the dump rope is unacceptable. This change induces additional vibrations in the rigging system and reduces the machine productivity.

Figure 4.2. Trajectories of hoist, dump, and drag ropes versus time during 1s of digging phase: (a) angular trajectories and (b) linear displacements of hoist and drag ropes

The integration of the DAEs in equation (4.2) with the proposed initial conditions,
(ππ_{ππ,1}) with ππ = 0, . . ,5, produced a singularity phenomenon shown in Figure 4.3. The second
initial values were chosen to minimize the residual values for integrating equation (4.1).
This apporach was further used in the calculations to improve the initial values of (ππ_{ππ,ππ}),
where n is the number of numerical experiments. At each experiment, the solver uses best

values to recalculate the residual at every time step and tries to mimimize it to zero. This kind of numerical study is relevant to an optimization method with inequality constraints.

Figure 4.3. Singularity of the kinematics model at time 1.52 s of digging phase

Another method for integrating the DAEs was to start with a zero initial value to every trajectory variable. The alogrithm with a predefined residual minimization embedded option tries to find better values and proceed to the final step. The output of the numerical experiment is shown in Figure 4.4 (a) and another singular location was found at time 7.16 seconds due to the stiffness of equations. Improvements in the numerical integration at time 7.16 seconds of 10-second total digging time do not permit the selection of initial values due to erroneous results in the linear displacements as seen in Figure 4.4 (b). Moreover, the field expirements show that a rapid increase in the trajectory of the hoist rope is not likely to happen during the digging phase. Thus, the numerical stablity of the solution algorithm is not beneficial and cannot be further expanded for finding the actual forces and torques in the front-end assembly.

**β’ Minimal set of constraint equations. The numerical analyses that were done a**
full set of constraint equations (4.2) resulted in singularity and untable solution process.
Numerical stability issues in these analyses are attributed to the fact that the functions of
linear displacements cause the integration to fail. In other words, those equations create
redundancy that can be avoided by eliminating them. The redundancy means that the
excavator has more DOFs available than the number required to do the task. This reason
has led to a reduction in the number of the constraint equations to a minimal set as shown
in equation (4.4). The numerical solution of these equations is also challenging due to the
presence of embedded constraints and the nonlinearities that follow the differentiations.
Equation (4.4) forms a set of DAEs that are used to find the initial conditions using
FindRoot or NDSolve, which are built-in subroutines in Mathematica. The same analogy
is applied to find the initial conditions of all kinematics quantities (ππ_{4}, ππ_{5}, ππ_{6}). It is
important to note that the differential kinematics constraint equations (4.4) may yield
multiple solutions or may not have any solutions depending the nature of the geometric
constraints.

Figure 4.4. Singularity of the kinematics model at time 7.16 s of digging phase: (a) angular trajectories and (b) linear displacements of hoist and drag ropes

1.32 c6+ 2.54 s4+ (75 + 2.54 t)c4qΜ4 β 1.715 s4 qΜ4+ 10.5 s5 qΜ5+
(75 β 1.32t)s_{6}qΜ_{6} = 0
2.54 c_{4}β 1.32 s_{6}β 1.715 c_{4} qΜ_{4}β (75 + 2.54 t)s_{4}qΜ_{4}β 10.5 c_{5}qΜ_{5}+
(75 β 1.32t)c6qΜ6 = 0
2.64 (75 β 1.32t) + 5.08 (75 + 2.54t) β53.34(75+2.54t) c[1.57+q_{οΏ½2.941β+(75+2.54t)}_{2}4+q5] β
290.4 c

### οΏ½

37Ο_{180}+ q

_{6}

### οΏ½

+21 οΏ½2.94β + (75 + 2.54t)2_{ s[1.57 + q}

_{4}

_{+ q}

_{5}

_{](qΜ}

_{4}

_{+ qΜ}

_{5}

_{) }β220(75 β 1.32t) s οΏ½37Ο

_{180}+ q6οΏ½ qΜ6 = 0

Equation (4.4) can be rewritten in a matrix form as shown in equations (4.5) and will have a solution only when the Jacobian matrix οΏ½J =ππFππ(qi)

ππqi οΏ½ is a full rank matrix.

οΏ½β1.715 c4(75 + 2.54 t)c4 β (75 + 2.54 t)s4β 1.715 s4 β10.5 c510.5 s5 (75 β 1.32t)s6(75 β 1.32t)c6 Z6 Z6 Z7 οΏ½ οΏ½qΜqΜ45 qΜ6 οΏ½ = οΏ½β1.32 cβ2.54 c46β 2.54 s+ 1.32 s64 Z8 οΏ½ (4.5)

The Jacobian is a 3x3 matrix given in the equation (4.6). Thus, the rank of the Jacobian
must be 3.
J = οΏ½β1.715 c(75 + 2.54 t)c4 β (75 + 2.54 t)s4β 1.715 s44 β10.5 c10.5 s55 (75 β 1.32t)c(75 β 1.32t)s66
Z6 Z6 Z7
οΏ½ (4.6)
Z6 = 21 οΏ½2.94β + (75 + 2.54t)2 s[1.57 + q4+ q5] (4.7)
Z7 = 220(75 β 1.32t) s οΏ½37Ο_{180}+ q6οΏ½ (4.8)
Z8= 2.64 (75 β 1.32t) + 5.08 (75 + 2.54t) β53.34(75+2.54t) c[1.57+q_{οΏ½2.941β+(75+2.54t)}_{2}4+q5]

### β290.4 c οΏ½37Ο 180+ q6οΏ½ (4.9) (4.4)

Figure 4.5. Trajectories of hoist, dump, and drag ropes versus time during digging phase

Integration of equation (4.4) was successfully performed using zero initial values
of the trajectory functions q_{4}, q_{5}and q_{6}, in Figure 4.5. The singularity behavior was not
an issue and the integration started and ended at the required interval of time [0, 10]
seconds. However, the evolution of the trajectory functions through this interval does not
meet the limits of the machineβs operational space. These responses are due to the fact that
the initial values used by the solver at time π‘π‘ = 0 second are still out of range. It can be
concluded that there is a trade-off between the initial value search and numerical stability
of the applied method. To circumvent this problem, subsection (4.1.1.3) highlights the use
of Baumgarte Stabilization Technique (BST) to solve a stiff DAE with inequality
constraints.