Global Energy-Optimal Redundancy Resolution of Hydraulic Manipulators: Experimental Results for a Forestry Manipulator

Abstract

This paper addresses the energy-inefficiency problem of four-degrees-of-freedom (4-DOF) hydraulic manipulators through redundancy resolution in robotic closed-loop controlled applications. Because conventional methods typically are local and have poor performance for resolving redundancy with respect to minimum hydraulic energy consumption, global energy-optimal redundancy resolution is proposed at the valve-controlled actuator and hydraulic power system interaction level. The energy consumption of the widely popular valve-controlled load-sensing (LS) and constant-pressure (CP) systems is effectively minimised through cost functions formulated in a discrete-time dynamic programming (DP) approach with minimum state representation. A prescribed end-effector path and important actuator constraints at the position, velocity and acceleration levels are also satisfied in the solution. Extensive field experiments performed on a forestry hydraulic manipulator demonstrate the performance of the proposed solution. Approximately 15–30% greater hydraulic energy consumption was observed with the conventional methods in the LS and CP systems. These results encourage energy-optimal redundancy resolution in future robotic applications of hydraulic manipulators.

1. Introduction

Hydraulic loaders are powerful, typically 4-degrees-of-freedom (4-DOF) manipulators, that are used for heavy-duty material handling on mobile machines, such as trucks and forest forwarders. Currently, hydraulic loaders are predominantly operated by humans, but loader manufacturers are interested in diversifying their offerings through robotic functionality to remain competitive. On the research side, and in the industry, subjects dealing with reducing reliance on the expertise of human operators in loader control have consequently received much attention. Several next-generation robotic semi-automated approaches, for example, have been proposed to help alleviate the operator’s workload by automating the repetitive work tasks that would normally be the operator’s responsibility [1,2]. On one hand, these semi-automated approaches discuss closed-loop control, and on the other hand, they involve redundancy resolution because hydraulic loaders are typically kinematically redundant. Although both problems are important, it is safe to estimate that redundancy resolution of hydraulic manipulators has received less and too little interest in academia in relation to the importance of redundancy resolution compared to closed-loop control problems. This is probably because so many different types of redundancy resolutions [3,4,5,6,7,8,9] have been proposed for general systems that the general understanding is that these redundancy resolutions can effectively reduce the energy consumption of hydraulic manipulators, for example, through the minimum velocity norm or weighted pseudo-inverse solutions [10]. One of the aims in this paper is to show that these general redundancy resolutions can be sub-optimal at the valve-controlled hydraulic system level.

Moreover, because there are many different redundancy resolution strategies for different objectives, choosing the most compelling objective for a hydraulic manipulator is not always straightforward. Looking at the handful of redundancy resolutions dealing with specific problems and objectives pertaining to hydraulic manipulators, however, it becomes evident that global energy-optimal redundancy resolutions are lacking. To explain the importance of global energy-optimal redundancy resolution, system energy inefficiency is a major problem in hydraulic manipulators [11]. Therefore, to address this problem, we resolve the redundancy in a highly desirable hydraulic energy-optimal manner and effectively treat the inherently limited joint motions of hydraulic manipulators by enforcing important actuator position, velocity and acceleration limits in our solution.

1.1. Literature Review

In the literature, Mettin et al. [12], Löfgren [13] and Ortiz Morales et al. [1] addressed the joint-limited minimum-time redundancy resolution of hydraulic manipulators, and Löfgren [13] has also studied redundancy resolution which maximises load capacity. Furthermore, Beiner and Mattila [10] have derived a pseudo-inverse method for redundancy resolution in the hydraulic manipulator’s actuator space to minimise the cylinder velocity norm or the weighted kinetic energy instantaneously. In addition, Flacco et al. [9] have proposed a joint-limited pseudo-inverse-based algorithm, which is also suitable for hydraulic applications, that saturates joint velocities in the null space. These different solutions, however, do not discuss hydraulic energy-optimality. In addition, generally minimising the energy consumption of actuators, as proposed by Vukobratovic and Kircanski [14] as a local solution and developed by Hirakawa and Kawamura [15] as an optimal control problem, can be a sub-optimal strategy for hydraulic manipulators. This is a direct consequence of the fact that the actuators’ energy consumption does not equate to the energy consumption of the typical valve-controlled load-sensing (LS) and constant-pressure (CP) hydraulic systems because the pressure losses encountered in these systems are not considered. It is important to derive energy-optimal redundancy resolutions for these commonly used systems. Whereas LS systems are predominantly used in hydraulic commercial manipulators in general, CP systems are used mostly in robotic closed-loop controlled applications and next-generation energy-efficient secondary-controlled applications because of the systems’ superior damping and simplicity compared with the LS systems [16,17,18].

1.2. Paper Contribution and Organization

In this paper, we minimise the energy consumption of the widely used LS and CP hydraulic systems through redundancy resolution to address the major problem of energy-inefficiency in hydraulic manipulators. This solution is in stark contrast to known conventional methods that are not energy-optimal at the valve-controlled actuator and hydraulic power system interaction level. To summarise, our redundancy resolution entails moving the manipulator hydraulic cylinders in an energy-optimal manner while the manipulator’s end-effector follows a prescribed path, and actuator limits at the position, velocity and acceleration levels are satisfied. To address both of these hydraulic systems optimally, we propose a tailored cost function of the energy consumption for each. For example, the simpler CP system is shown to require minimisation of the pump flow rate. This objective leads to our effective solution in which the cylinder’s effective areas scale the cylinder velocities in a discontinuous manner. Minimising the CP system energy cost, hence, intuitively amounts to minimising the cylinder volumes displaced in the solution. The cylinder areas among the manipulator’s hydraulic cylinders, including the different piston and rod-side areas of the asymmetric cylinders, render minimisation of the CP system’s energy consumption, in particular, a highly tractable objective.

With the objectives in place, the redundancy could obviously be resolved by formulating the optimal control problem so that the motions of all the joints are optimised subject to specific motion constraints. However, that approach would result in a high-dimensional problem, which would be unsuitable for a global dynamic programming (DP) solution and generally inconvenient. To reduce the problem dimensionality and computational complexity, in this study, only the motion of the redundant extension cylinder is directly optimised. The optimal non-redundant joint motions are, instead, resolved by employing non-singular inverse kinematics, motivated by the formulation used by Choi [19] and Flacco et al. [20]. Moreover, in the present study, hydraulic system properties are incorporated in the cost functions proposed in steady-state equations. This approach results in our desirable minimum-state problem, which requires only two states and one control and which uses the computational resources available much more sparingly than the approach of simultaneously optimising all of the joint motions. However, because of the modelling simplifications in our problem formulation, a natural concern that could arise is whether significant hydraulic aspects of the problem are neglected. Thus, extensive experiments, which compare the energy consumptions of the redundancy resolutions proposed in this study against other global and instantaneous strategies in the LS and CP systems is presented.

This energy-optimal solution is primarily suitable as an offline approach, much like the minimum-time solutions that were proposed by Ortiz Morales et al. [1] and Mettin et al. [12], and it is practical as such. Extensive measurements by Ortiz Morales et al. [1], Mettin et al. [12] and Westerberg [2] have shown that the typical workspace activity for 4-DOF manipulator is concentrated on specific areas, which consequently means that there are in practice a limited number of significant point-to-point movements to optimise with respect to energy consumption. However, at the implementation stage, it might be feasible to pre-compute the energy-optimal trajectories over a few time horizons for all of these significant point-to-point movements so that it would be possible to choose between slower and faster energy-optimal trajectories in real-time. Furthermore, with our proposed methodology, one significant implication is that only the redundant joint’s trajectories are necessary for real-time applications. If having more storage space were considered more important than having more computational power reserve, then the non-redundant motion trajectories could be resolved online from the remaining non-singular inverse kinematics.

This paper is an extension of our previous work [21]. As an extension of previous work, redundancy resolution is expanded from vertical plane three-degrees-of-freedom (3-DOF) serial-link hydraulic manipulators to complex 4-DOF serial-link hydraulic manipulators, which include boom base motion. The system model proposed for the energy-optimal solution is also systematically presented in this study. A third contribution of this paper compared to our previous simulation study concerning 3-DOF serial-link hydraulic manipulators is the extensive field experimental evaluation of the redundancy resolution proposed for 4-DOF serial-link hydraulic manipulators and its comparison to existing conventional methods. In this experimental evaluation, the performance of the proposed DP solution is compared to a “bootstrapped” DIDO solution [22] by employing the proposed problem formulation. The performance of the globally hydraulic energy-optimised redundancy resolution has not been previously experimentally evaluated.

This paper is organized as follows. First, the test system, which is a commercial forestry hydraulic manipulator, is described in Section 2. Second, a simplified mathematical model of this complex manipulator is presented in Section 3. In Section 4, the hydraulic manipulator’s optimal control problem is reduced to the bare essentials with the minimum-state formulation. The problem’s desirable global dynamic programming solution is discussed in Section 5. Finally, we experimentally evaluate the proposed redundancy resolution by comparing its performance to that of conventional methods in Section 6. The experimental results are from the real test-bed introduced in Section 2. Conclusions are provided in Section 7.

2. Test-Bed

Our test-bed is the K70M hydraulic manipulator on the Ponsse Caribou S10 8-wheel forest forwarder (Figure 1). In the following, we describe this manipulator in more detail.

2.1. System Description

The K70M hydraulic manipulator is a standard log-loader that has a grapple end-effector attached to the crane tip. The manipulator has four principal joints, which consist of the base, lift, tilt and extension joints shown in Figure 1. These joints are actuated by using proportional valve-controlled hydraulic cylinders (see

Table 1. Hydraulic cylinder dimensions of the test-bed manipulator.

i	$\oslash {\mathit{A}}_{\mathbf{A}\mathit{i}}/{\mathit{A}}_{\mathbf{B}\mathit{i}}$ − Stroke [mm]	Explanation
1	$2\times \oslash 115/115-720$	Base cylinders
2	$\oslash 120/70-610$	Lift cylinder
3	$\oslash 115/60-725$	Tilt cylinder
4	$\oslash 70/45-1850$	Extension cylinder

). The manipulator also has two telescoping links in series to extend its reachability. Measuring from the base to the tip, the manipulator can extend to approximately 9.1 m, and at this distance, its load-carrying capacity is 650 kg.

A two-stage, pressure-compensated Parker K170LS mobile hydraulic valve is used to control the manipulator’s cylinders; see Figure 2, where the valve’s spool types are illustrated. The mobile hydraulic valve also includes very complex in-built hydro-mechanical functionality, such as a standard counter pressure valve, which raises the tank pressure, and anti-cavitation valves in each actuator port (not shown in Figure 2). In the over-running load case, in particular, these anti-cavitation valves allow any fluid flow rate from the tank line to the non-loaded cylinder chamber to exist before the chamber pressure drops to a level below a preset minimum hydraulic cylinder pressure. This hydro-mechanical functionality can be significant since these anti-cavitation valves may reduce pump flow rate requirement.

The manipulator, as an industrial standard, is powered by a LS hydraulic pump; see [23,24,25]. Therefore, an electro-hydraulic CP system was built in parallel to allow comparison to the LS system for the experiments. Hence, the same hydraulic pump and pump LS regulator valve can be used to control the pump displacement. The hydraulic pump is an A11VLO130DRS variable axial piston pump from Bosch Rexroth with a maximum flow rate of 325 $\mathrm{d}{\mathrm{m}}^3/\min$ at 2500 rotations per minute (including the charge pump flow) and 21.5 MPa maximum supply pressure.

The LS system’s conventional low-level operating principle, which is reported in [23], is illustrated in Figure 3a. In this system, the maximum driven chamber pressure among the base, lift, tilt and extension cylinders is fed back to the variable displacement pump’s hydro-mechanical controller, which sets the pump supply pressure a LS pressure margin $\Delta p_{\mathrm{LS}}$ above this feedback value. With the CP electro-hydraulic control method of the hydro-mechanical feedback system shown in Figure 3b, on the other hand, the pressure reference for the variable displacement pump’s hydro-mechanical controller is generated directly from the supply pressure based on the electrical input of the proportional pressure relief valve between the minor (0.7 mm) orifice and the tank-connected pressure valve. This control method is called “remote control” in the Bosch Rexroth catalog [26]. We have simplified the symbol of the metre-in side load control valve and omitted the metre-out side for clarity.

2.2. Instrumentation and Data Acquisition

The hydraulic manipulator was equipped with multiple sensors for the experiments:

(1)

the base angle was measured using an RCC absolute angle resolver from LTN with a resolution of 0.01${}^{\circ }$;

(2)

the lift and tilt angles were measured using incremental angle encoders ROD 456–5000 from Heidenhain with resolutions of 0.0072${}^{\circ }$ and 0.0029${}^{\circ }$, after 10-fold and 25-fold interpolation;

(3)

the extension cylinder position was measured with a DG 60 L incremental position encoder from Stegmann with a resolution of 160 $\mathsf{\mu }$m;

(4)

the pump supply pressure and, for identification purposes, cylinder pressures were measured using an NAH250 analogue pressure transducer from Trafag with an accuracy of ±0.1% from full-scale 25 MPa and

(5)

the pump flow rate was measured using a gear type flow meter VC 5 F2 PS from Kracht with a resolution of 0.005 $\mathrm{d}{\mathrm{m}}^3$ and a range of −250 to 250 $\mathrm{d}{\mathrm{m}}^3$ per minute.

The joint encoder installations are shown in Figure 4. These sensor measurements were recorded at a 2 ms sampling rate, and the hydraulic valves were controlled via CAN using a dSPACE MicroAutoBox II processing unit for the experimental evaluation.

3. Model of the Kinematically Redundant Four-Degrees-of-Freedom Hydraulic Manipulator

Next, the test-bed is described in more detail. For redundancy resolution of this test-bed manipulator, the mathematical model of the manipulator is systematically presented in the following, with particular attention paid to the applicability of the model for optimising energy consumption.

3.1. Forward Kinematics in the Joint Space

The hydraulic manipulator is shown in Figure 5. The joint-position vector of this manipulator can be defined as follows:

$$\mathbf{q}={\left[\begin{array}{cccc}{q}_1& q_2& q_3& q_4\end{array}\right]}^{\mathrm{T}}$$

(1)

where $q_1$ denotes the base angle, $q_2$ denotes the lift angle, $q_3$ denotes the tilt angle and $q_4$ is the total extension length (a prismatic joint). These are the principal hydraulically actuated joints. All $q_i$ are taken positive counter-clockwise.

The manipulator’s end-effector position can be expressed from the base frame 0 to frame 4 in accordance with the classical Denavit-Hartenberg (DH) convention [27]. By using the DH transformation matrix ${\mathbf{A}}_i^{i-1}$, which determines the coordinate transformation between successive frames, we get:

$${\mathbf{A}}_4^0={\mathbf{A}}_1^0{\mathbf{A}}_2^1{\mathbf{A}}_3^2{\mathbf{A}}_4^3$$

(2)

as the overall coordinate transformation expressed in the base frame. The transformation matrix ${\mathbf{A}}_i^{i-1}\in {\mathbb{R}}^{4x4}$ used consecutively is defined as:

$${\mathbf{A}}_i^{i-1}=\left[\begin{array}{cccc}{\mathrm{c}}_{{\theta }_i}& -{\mathrm{s}}_{{\theta }_i}{\mathrm{c}}_{{\alpha }_i}& {\mathrm{s}}_{{\theta }_i}{\mathrm{s}}_{{\alpha }_i}& a_i{\mathrm{c}}_{{\theta }_i}\\ {\mathrm{s}}_{{\theta }_i}& {\mathrm{c}}_{{\theta }_i}{\mathrm{c}}_{{\alpha }_i}& -{\mathrm{c}}_{{\theta }_i}{\mathrm{s}}_{{\alpha }_i}& a_i{\mathrm{s}}_{{\theta }_i}\\ 0& {\mathrm{s}}_{{\alpha }_i}& {\mathrm{c}}_{{\alpha }_i}& d_i\\ 0& 0& 0& 1\end{array}\right]$$

(3)

where, for example, ${\mathrm{c}}_{{\alpha }_i}$ denotes $cos\left({\alpha }_i\right)$ and ${\mathrm{s}}_{{\alpha }_i}$ denotes $sin\left({\alpha }_i\right)$. The DH parameters given in

Table 2. Denavit-Hartenberg parameters of the 4-DOF hydraulic manipulator.

Frame No.	${\mathit{a}}_{\mathit{i}}\left[\mathbf{m}\right]$	${\mathit{\alpha }}_{\mathit{i}}\left[\mathbf{rad}\right]$	${\mathit{d}}_{\mathit{i}}\left[\mathbf{m}\right]$	${\mathit{\theta }}_{\mathit{i}}\left[\mathbf{rad}\right]$
1	$o_{\mathrm{x}}=-0.208$	$\pi /2$	$o_{\mathrm{y}}=1.605$	$q_1$
2	$L_1=3.504$	0	0	$q_2$
3	$L_2=0.376$	$\pi /2$	0	$q_3+\pi /2$
4	0	$-\pi /2$	$L_3+q_4=2.13+q_4$	0

determine these matrix elements.

By transforming from the base frame to the $x_{\mathrm{w}}$, $y_{\mathrm{w}}$ and $z_{\mathrm{w}}$ coordinates, which denote the x-, y- and z-positions of the end-effector expressed in the global frame, respectively, we get:

$$x_{\mathrm{w}}=cos\left(q_1\right)\left\{L_{1}cos\left(q_2\right)-L_{2}sin(q_2+q_3)+(L_3+q_4)cos(q_2+q_3)+o_{\mathrm{x}}\right\}$$

(4)

$$y_{\mathrm{w}}=L_{1}sin\left(q_2\right)+L_{2}cos(q_2+q_3)+(L_3+q_4)sin(q_2+q_3)+o_{\mathrm{y}}\mathrm{and}$$

(5)

$$z_{\mathrm{w}}=sin\left(q_1\right)\left\{L_{1}cos\left(q_2\right)-L_{2}sin(q_2+q_3)+(L_3+q_4)cos(q_2+q_3)+o_{\mathrm{x}}\right\}$$

(6)

as the final end-effector position, which we denote with ${\mathbf{x}}_{\mathrm{t}}={\begin{array}{c}\left[x_{\mathrm{w}}{y}_{\mathrm{w}}{z}_{\mathrm{w}}\right]\end{array}}^{\mathrm{T}}$. Constants $L_1$ and $L_3$ denote the manipulator link lengths, $L_2$ is an offset between the tilt joint and the extension link and $o_{\mathrm{x}}$ and $o_{\mathrm{y}}$ are the offsets in the direction of the x- and y-axes, respectively, from the global frame to the lift joint. The workspace of the test-bed manipulator in the vertical plane is depicted in Figure 6 by using the parameters given in Table 2.

Differentiating ${\mathbf{x}}_{\mathrm{t}}$ in view of Equations (4)–(6) with respect to time yields the end-effector velocity:

$${\dot{\mathbf{x}}}_{\mathrm{t}}=\mathbf{J}\left(\mathbf{q}\right)\dot{\mathbf{q}}={\mathbf{J}}_1{\dot{q}}_1+{\mathbf{J}}_2{\dot{q}}_2+{\mathbf{J}}_3{\dot{q}}_3+{\mathbf{J}}_4{\dot{q}}_4$$

(7)

where ${\dot{\mathbf{x}}}_{\mathrm{t}}={\begin{array}{c}\left[{\dot{x}}_{\mathrm{w}}{\dot{y}}_{\mathrm{w}}{\dot{z}}_{\mathrm{w}}\right]\end{array}}^{\mathrm{T}}$, the joint velocity vector $\dot{\mathbf{q}}={\begin{array}{c}\left[{\dot{q}}_1{\dot{q}}_2{\dot{q}}_3{\dot{q}}_4\right]\end{array}}^{\mathrm{T}}$ and the Jacobian matrix $\mathbf{J}\left(\mathbf{q}\right)$ is the partial derivative of ${\mathbf{x}}_{\mathrm{t}}$ with respect to $\mathbf{q}$. This equation was also expressed by using the column vectors ${\mathbf{J}}_i$, for $i\in \{1,2,3,4\}$, for future reference. Differentiating ${\dot{\mathbf{x}}}_{\mathrm{t}}$ once more yields:

$${\ddot{\mathbf{x}}}_{\mathrm{t}}=\dot{\mathbf{J}}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}+\mathbf{J}\left(\mathbf{q}\right)\ddot{\mathbf{q}}$$

(8)

with the end-effector acceleration ${\ddot{\mathbf{x}}}_{\mathrm{t}}$ and the joint acceleration vector $\ddot{\mathbf{q}}$ defined similarly to Equation (7). The elements of the Jacobian time derivative $\dot{\mathbf{J}}$ can be resolved from:

$$\dot{J}{(\mathbf{q},\dot{\mathbf{q}})}_{ji}=\frac{\partial J{\left(\mathbf{q}\right)}_{ji}}{\partial \mathbf{q}}\dot{\mathbf{q}}=\sum _{n=1}^4\frac{\partial J{\left(\mathbf{q}\right)}_{ji}}{\partial q_n}{\dot{q}}_n$$

(9)

over the workspace dimensions $j\in \{1,2,3\}$ and the joint dimensions $i\in \{1,2,3,4\}$. Because the Jacobian matrix $\mathbf{J}\left(\mathbf{q}\right)\in {\mathbb{R}}^{3\times 4}$ is non-square, the inverse kinematics clearly have an infinite number of solutions inside the workspace. Thus, a redundancy resolution is required.

3.2. Joint Space to Actuator Space Transformations

Transformations, which yield the cylinder coordinates as a function of the joints, are presented next. By using the chain rule, we obtain:

$$c_i=c_i\left(q_i\right),$$

(10)

$$v_i=\frac{\partial c_i\left(q_i\right)}{\partial q_i}{\dot{q}}_i=r_i{\dot{q}}_i\mathrm{and}$$

(11)

$$a_i=\frac{\partial r_i}{\partial q_i}{\dot{q}}_i^2+r_i{\ddot{q}}_i$$

(12)

as the general equations, where $c_i\left(q_i\right)$, for the actuator index $i\in \{1,2,3,4\}$, is the cylinder length, $v_i$ is the cylinder (piston) velocity, $a_i$ is the cylinder (piston) acceleration, not to be confused with the DH parameter, and $r_i$ is the cylinder torque arm. The joint space to actuator space transformations of the test-bed are presented in the following.

Firstly, the base cylinder has a constant torque arm $r_1$ because a rack-and-pinion mechanism, as illustrated in Figure 7, is used (Equation (11)). Furthermore, the torque arm $r_2$ of the lift cylinder can be written as:

$$\begin{array}{cc}\hfill {\beta }_{13}& =q_2-{\beta }_{11}-{\beta }_{12}+\pi /2\hfill \\ \hfill r_2\left(q_2\right)& =\frac{L_{11}{L}_{12}sin\left({\beta }_{13}\right)}{c_2}=\frac{L_{11}{L}_{12}sin\left({\beta }_{13}\right)}{\sqrt{L_{11}^2+L_{12}^2-2L_{11}{L}_{12}cos\left({\beta }_{13}\right)}}\hfill \end{array}$$

(13)

as a function of the lift angle $q_2$ and the lift cylinder length $c_2$. The constants $L_{11}$, $L_{12}$, ${\beta }_{11}$ and ${\beta }_{12}$ and the joint angle $q_2$ are illustrated in Figure 8a.

The four-bar linkage of the tilt cylinder is governed by more complex kinematic equations. These equations, from the tilt angle $q_3$ to tilt cylinder length $c_3$, can be written as follows:

$$\begin{array}{cc}\hfill {\alpha }_{23}& =\pi +q_3+{\alpha }_{21}+{\alpha }_{22}\hfill \\ \hfill L_{26}& =\sqrt{L_{22}^2+L_{23}^2-2L_{22}{L}_{23}cos\left({\alpha }_{23}\right)}\hfill \\ \hfill {\beta }_{23}& =\mathrm{asin}\left(L_{23}/L_{26}sin\left({\alpha }_{23}\right)\right)\hfill \\ \hfill {\beta }_{24}& =\mathrm{acos}\left((L_{25}^2+L_{26}^2-L_{24}^2)/\left(2L_{25}{L}_{26}\right)\right)\hfill \\ \hfill {\beta }_{22}& =\pi -{\beta }_{21}-{\beta }_{24}+{\beta }_{23}+{\alpha }_{21}\hfill \\ \hfill c_3\left(q_3\right)& =\sqrt{L_{21}^2+L_{25}^2-2L_{21}{L}_{25}cos\left({\beta }_{22}\right)}\hfill \\ & \approx s_5{q}_3^5+s_4{q}_3^4+s_3{q}_3^3+s_2{q}_3^2+s_1{q}_3+s_0\hfill \\ & =f_q(\mathbf{s},q_3),\hfill \end{array}$$

(14)

where ${\alpha }_{21}$, ${\alpha }_{22}$, ${\beta }_{21}$, $L_{21}$, $L_{22}$, $L_{23}$, $L_{24}$ and $L_{25}$ are the constants and ${\beta }_{22}$, ${\beta }_{23}$, ${\beta }_{24}$ and $L_{26}$ are the variables (see Figure 8b, where $q_3$ is negative). Laws of cosines and sines were used. A quintic function $f_q(\mathbf{s},q_3)$, with $\mathbf{s}\in {\mathbb{R}}^6$ as the coefficient vector, can be fitted to the $c_3$ data to simplify the equation of the tilt torque arm $r_3$. Hence, differentiating $c_3$ with respect to time in view of Equation (11) yields:

$$\begin{array}{cc}\hfill v_3& =-\frac{\partial f_q(\mathbf{s},q_3)}{\partial q_3}{\dot{q}}_3\hfill \\ & =-(5s_5{q}_3^4+4s_4{q}_3^3+3s_3{q}_3^2+2s_2{q}_3+s_1){\dot{q}}_3\hfill \\ & =-r_3\left(q_3\right){\dot{q}}_3\hfill \end{array}$$

(15)

as the tilt cylinder’s piston velocity. Multiplication by minus one considers the opposite signs of the piston and angular velocities. The “torque arm” $r_4$ of the telescoping extension mechanism is $1/2$ because the extension cylinder velocity is half of the joint velocity due to the telescope’s mechanical multiplier.

3.3. Cylinder Forces

Manipulator dynamics are considered through the hydraulic cylinder forces, which can be resolved from the inverse rigid-body dynamics as follows:

$$\mathbf{F}={\mathbf{R}}^{-1}(\underset{\mathit{\tau }}{\underbrace{\mathbf{M}\left(\mathbf{q}\right)\ddot{\mathbf{q}}+\mathbf{G}\left(\mathbf{q}\right)+{\mathbf{F}}_{\mathrm{f}}(\mathbf{q},\dot{\mathbf{q}})}})$$

(16)

where the cylinder force vector $\mathbf{F}$ equals ${\begin{array}{c}\left[F_1{F}_2{F}_3{F}_4\right]\end{array}}^{\mathrm{T}}$, $\mathit{\tau }$ is the corresponding cylinder torque vector, $\mathbf{M}\left(\mathbf{q}\right)\in {\mathbb{R}}^{4x4}$ denotes the positive definite, symmetric inertia matrix at joint position $\mathbf{q}$, $\mathbf{G}\left(\mathbf{q}\right)\in {\mathbb{R}}^4$ denotes the gravitational torque vector and ${\mathbf{F}}_{\mathrm{f}}(\mathbf{q},\dot{\mathbf{q}})\in {\mathbb{R}}^4$ includes the hydraulic cylinder and joint friction forces. The inertia matrix and the gravitational load vector can be resolved via well-known Lagrangian conventions by summing the dynamical contributions of each link in the global frame [27]. Resolving these components regarding our test-bed is mostly standard practice, with the exception that the extension cylinder’s gravitational component must be doubled to consider the mechanical disadvantage of the telescope mechanism (The mechanism doubles the extension cylinder stroke). Because the centrifugal and Coriolis forces are in low magnitude, their effects can be omitted. Conventional exponential friction force models can be primarily exploited (see [28]), but the extension joint’s friction force is shown to be accurately represented by using the developed unconventional model (see Appendix B for the test-bed’s cylinder force modelling results).

The inverse torque arm matrix is defined as in:

$${\mathbf{R}}^{-1}=\left[\begin{array}{cccc}1/r_1& 0& 0& 0\\ 0& 1/r_2& 0& 0\\ 0& 0& 1/r_3& 0\\ 0& 0& 0& 1/r_4\end{array}\right]$$

(17)

Consequently, these positive matrix elements, denoted by $r_i^{-1}$, are the inverses of the hydraulic cylinder torque arms. However, $r_4^{-1}$ is a non-physical multiplier because the fourth joint is prismatic.

3.4. Steady-State Cylinder Pressures Using a Mobile Hydraulic Valve

Finally, moving on to the hydraulic system equations, we introduce simplifications to keep the system dimensionality at a manageable level for optimisation purposes. To introduce our subsequent simplifications, firstly, consider the pressure dynamics of the hydraulic cylinders denoted by [29]:

$${\dot{p}}_{\mathrm{A}}=\frac{B_{\mathrm{eff}}}{V_{0\mathrm{A}}+x_{\mathrm{p}}{A}_{\mathrm{A}}}\left(Q_{\mathrm{A}}-v{A}_{\mathrm{A}}\right)\hspace{0.17em}\mathrm{and}$$

(18)

$${\dot{p}}_{\mathrm{B}}=\frac{B_{\mathrm{eff}}}{V_{0\mathrm{B}}+(L_{\max }-x_{\mathrm{p}})A_{\mathrm{B}}}\left(Q_{\mathrm{B}}+v{A}_{\mathrm{B}}\right)$$

(19)

where $B_{\mathrm{eff}}$ is the effective bulk modulus, $V_{0\mathrm{A}}$ and $V_{0\mathrm{B}}$ are the dead volumes, $x_{\mathrm{p}}$ is the piston position, v is the piston velocity, $A_{\mathrm{A}}$ is the piston-side area, $A_{\mathrm{B}}$ is the rod-side area and $Q_{\mathrm{A}}$ and $Q_{\mathrm{B}}$ are the flow rates to the cylinder chambers. Actuator indices were omitted for brevity. Solving the steady-state cylinder velocity from Equations (18) and (19) yields:

$$\begin{array}{c}\hfill v_{\mathrm{ss}}=Q_{\mathrm{A}}/A_{\mathrm{A}}=-Q_{\mathrm{B}}/A_{\mathrm{B}}\end{array}$$

(20)

Now, we omit these hydraulic system dynamics and focus on deriving the hydraulic system’s steady-state pressure equations by utilising Equation (20). Because the flow rates to the hydraulic cylinders are controlled by using the mobile hydraulic valve’s two vastly different spools, as shown in Figure 2, we derive the steady-state cylinder pressures in both cases.

3.4.1. Typical Pressure-Compensated Non-Differential Spool Valve

Firstly, for the asymmetric spool valve shown in Figure 2a, which is another typical pressure-compensated, hydro-mechanical valve, the turbulent flow equations are given by [29]:

$$Q_{\mathrm{A}}=\left\{\begin{array}{cc}{K}_{\mathrm{vPA}}\left(I\right)\sqrt{\Delta p_{\mathrm{c}}},\hfill & \mathrm{if}\mathrm{I}\ge 0\left(\mathrm{cylinder}\mathrm{extends}\right)\hfill \\ -K_{\mathrm{vAT}}\left(I\right)\sqrt{p_{\mathrm{A}}-p_{\mathrm{T}}},\hfill & \mathrm{if}\mathrm{I}<0\left(\mathrm{cylinder}\mathrm{retracts}\right)\hfill \end{array}\right.$$

(21)

$$Q_{\mathrm{B}}=\left\{\begin{array}{cc}-K_{\mathrm{vBT}}\left(I\right)\sqrt{p_{\mathrm{B}}-p_{\mathrm{T}}},\hfill & \mathrm{if}\mathrm{I}\ge 0(\mathrm{cyl}.\mathrm{ext}.)\hfill \\ K_{\mathrm{vPB}}\left(I\right)\sqrt{\Delta p_{\mathrm{c}}},\hfill & \mathrm{if}\mathrm{I}<0(\mathrm{cyl}.\mathrm{ret}.)\hfill \end{array}\right.$$

(22)

where I is the valve’s control current, $\Delta p_c$ is the constant pressure drop across the main spool due to the pressure compensation, $p_{\mathrm{T}}$ is the tank pressure and $K_{\mathrm{vPA}}$, $K_{\mathrm{vAT}}$, $K_{\mathrm{vBT}}$ and $K_{\mathrm{vPB}}$ are the flow coefficients of the P-A, A-T, B-T and P-B control edges, respectively. The flow coefficients $K_{\mathrm{vPA}}\left(I\right)$ and $K_{\mathrm{vBT}}\left(I\right)$ are monotonically increasing positive functions for a valve current I greater than a positive dead zone current $I_{\mathrm{DZ}}^{(+)}$. Otherwise, they have a value of zero. The flow coefficients $K_{\mathrm{vAT}}\left(I\right)$ and $K_{\mathrm{vPB}}\left(I\right)$ are similarly defined as monotonically increasing positive functions for a current I less than a negative dead zone current $I_{\mathrm{DZ}}^{(-)}$. Here, ideal pressure compensation was assumed, and reversed flow was neglected to simplify analysis.

Utilising Equations (20)–(22) and the external force (which equals the cylinder force $F=p_{\mathrm{ssA}}{A}_{\mathrm{A}}-p_{\mathrm{ssB}}{A}_{\mathrm{B}}$) to resolve the steady-state pressures yields [30]:

$$\begin{array}{c}\hfill p_{\mathrm{ssA}}=\left\{\begin{array}{cc}F/A_{\mathrm{A}}+\frac{A_{\mathrm{B}}}{A_{\mathrm{A}}}{p}_{\mathrm{ssB}},\hfill & \mathrm{for}\mathrm{I}>I_{\mathrm{DZ}}^{(+)}\hfill \\ \underset{\mathrm{Coeff}.}{\underbrace{{\left(\frac{K_{\mathrm{vPB}}{A}_{\mathrm{A}}}{K_{\mathrm{vAT}}{A}_{\mathrm{B}}}\right)}^2}}\Delta p_{\mathrm{c}}+p_{\mathrm{T}},\hfill & \mathrm{for}\mathrm{I}<I_{\mathrm{DZ}}^{(-)}\hfill \end{array}\right.\end{array}$$

(23)

$$\begin{array}{c}\hfill p_{\mathrm{ssB}}=\left\{\begin{array}{cc}\stackrel{\mathrm{Coeff}.}{\overbrace{{\left(\frac{K_{\mathrm{vPA}}{A}_{\mathrm{B}}}{K_{\mathrm{vBT}}{A}_{\mathrm{A}}}\right)}^2}}\Delta p_{\mathrm{c}}+p_{\mathrm{T}},\hfill & \mathrm{for}\mathrm{I}>I_{\mathrm{DZ}}^{(+)}\hfill \\ -F/A_{\mathrm{B}}+\frac{A_{\mathrm{A}}}{A_{\mathrm{B}}}{p}_{\mathrm{ssA}},\hfill & \mathrm{for}\mathrm{I}<I_{\mathrm{DZ}}^{(-)}\hfill \end{array}\right.\end{array}$$

(24)

Interestingly, the cylinder back-pressures can be approximated by multiplying the pressure compensator’s constant pressure drop with a momentarily constant coefficient (labelled). These equations apply regardless of the load force direction. If the cylinder back-pressure is non-existent, as in single-acting lift cylinders, then these steady-state pressure approximations can be easily simplified.

3.4.2. Typical Pressure-Compensated Differential Spool Valve

The second spool type is the pressure-compensated differential spool valve (Figure 2b). Because of the spool’s differential position, the fluid flows are regenerative from the cylinder chamber B to cylinder chamber A via the valve’s P port in the cylinder extension case. This differential position reduces $Q_{\mathrm{p}}$, the required pump flow rate, to $Q_{\mathrm{A}}-\left|Q_{\mathrm{B}}\right|$. Hence, the actuator flow rates can be written as:

$$Q_{\mathrm{A}}=\underset{Q_{\mathrm{p}}}{\underbrace{K_{\mathrm{vPA}}\left(I\right)\sqrt{\Delta p_{\mathrm{c}}}}}+\underset{|Q_{\mathrm{B}}|}{\underbrace{|K_{\mathrm{vBP}}\left(I\right)|\sqrt{p_{\mathrm{B}}-p_{\mathrm{A}}-\Delta p_{\mathrm{c}}}}}$$

(25)

$$Q_{\mathrm{B}}=-K_{\mathrm{vBP}}\left(I\right)\sqrt{p_{\mathrm{B}}-p_{\mathrm{A}}-\Delta p_{\mathrm{c}}}$$

(26)

where $\mathrm{I}>0$ and $K_{\mathrm{vPA}}$ and $K_{\mathrm{vBP}}$ are defined as positive functions of current, $\Delta p_{\mathrm{c}}$ is a constant-pressure drop over the control edge P-A due to the pressure compensation and the sum of $p_{\mathrm{A}}$ and $\Delta p_{\mathrm{c}}$ equals the pressure in the valve’s P port, omitting some pressure losses ($Q_{\mathrm{B}}$ denotes the flow rate over the control edge B-P). Using Equations (20), (25) and (26), the steady-state cylinder chamber pressures can be resolved, yielding:

$$\begin{array}{c}\hfill p_{\mathrm{ssB}}=\frac{F}{A_{\mathrm{A}}-A_{\mathrm{B}}}+\left(\frac{K_{\mathrm{vPA}}^2{A}_{\mathrm{B}}^2{A}_{\mathrm{A}}}{K_{\mathrm{vPB}}^2{(A_{\mathrm{A}}-A_{\mathrm{B}})}^3}+\frac{A_{\mathrm{A}}}{A_{\mathrm{A}}-A_{\mathrm{B}}}\right)\Delta p_{\mathrm{c}}\end{array}$$

(27)

$$\begin{array}{c}\hfill p_{\mathrm{ssA}}=F/A_{\mathrm{A}}+p_{\mathrm{ssB}}\frac{A_{\mathrm{B}}}{A_{\mathrm{A}}}\end{array}$$

(28)

Here, the flow coefficients and the external force are again momentarily constant. Pressure losses in pipes, hoses and filters were omitted. In contrast to the non-differential case, the cylinder back-pressure $p_{\mathrm{ssB}}$ depends on the cylinder force. None of these derived steady-state pressures, however, depend explicitly on the supply pressure because of the pressure compensation.

4. Problem Formulation

Let us formulate our inherently continuous-time optimal control problem in discrete-time for numerical solution. By dividing the fixed optimisation time $t_{\mathrm{f}}$, which corresponds to the desired duration of the end-effector path, into N intervals of length $t_{\mathrm{f}}/N$, our optimal control problem can be expressed in discrete-time as follows:

$$\underset{{\pi }_u\in {\pi }_{\mathbb{U}}}{min}\sum _{k=0}^{N-1}{P}_{\mathrm{h},k}({\mathbf{q}}_k,{\dot{\mathbf{q}}}_k,{\ddot{\mathbf{q}}}_k)+{\ell }_k({\mathbf{q}}_k,{\dot{\mathbf{q}}}_k,{\ddot{\mathbf{q}}}_k)$$

(29)

subject to:

$$\begin{array}{cc}\hfill {\overline{\mathbf{x}}}_{k+1}={\overline{\mathbf{F}}}_k& ={\overline{\mathbf{x}}}_k+T_{\mathrm{s}},{\mathbf{F}}_k({\overline{\mathbf{x}}}_k,{\overline{u}}_k)\hfill \end{array}$$

(30)

where ${\pi }_u$ denotes the control policy and ${\pi }_{\mathbb{U}}$ denotes the admissible control policies. The hydraulic power $P_{\mathrm{h},k}$, defined at each time stage $k\in \{0,1,2,\dots ,N-1\}$, is the principal cost to be minimised, and the term ${\ell }_k$, which can be replaced by standard inequality constraints, penalises violations of actuator limits. For the performance cost $P_{\mathrm{h},k}$, effective functions which minimise hydraulic energy consumption are proposed because the conventional cost functions may be sub-optimal.

In the above, the system dynamics have been discretised using the Euler method and $T_{\mathrm{s}}$ as the integration step corresponding to the interval length $t_{\mathrm{f}}/N$. The state vector ${\overline{\mathbf{x}}}_k={\begin{array}{c}\left[{\overline{c}}_{4,k}{\overline{v}}_{4,k}\right]\end{array}}^{\mathrm{T}}$, with initial state ${\overline{\mathbf{x}}}_0$, contains the extension cylinder position and velocity. The extension cylinder acceleration ${\overline{a}}_{4,k}$ is used as the control ${\overline{u}}_k$. To optimise movements of the extension cylinder at the acceleration level, we create a double integrator system ${\mathbf{F}}_k({\overline{\mathbf{x}}}_k,{\overline{u}}_k)$, denoted with:

$$\begin{array}{c}\hfill {\mathbf{F}}_k({\overline{\mathbf{x}}}_k,{\overline{u}}_k)=\left[\begin{array}{cc}0& v_{4,\max }/c_{4,\max }\\ 0& 0\end{array}\right]{\overline{\mathbf{x}}}_k+\left[\begin{array}{c}0\\ a_{4,\max }/v_{4,\max }\end{array}\right]{\overline{u}}_k\end{array}$$

(31)

where the maximum extension cylinder position, velocity and acceleration, $c_{4,\max }$, $v_{4,\max }$ and $a_{4,\max }$, respectively, have been applied to normalise the physical units for the cylinder’s position, velocity and acceleration. These normalised units can be transformed into conventional physical units by employing scaling equations ${\mathbf{x}}_k=\mathrm{diag}(c_{4,\max },v_{4,\max }){\overline{\mathbf{x}}}_k$ and $u_k=a_{4,\max }{\overline{u}}_k$. These extension cylinder coordinates can also be converted into joint coordinates by multiplying the cylinder motion by $1/r_4$, in view of Equations (10)–(12). This final transformation is required for the next stage.

4.1. Using Inverse Kinematics to Resolve the Motion of Non-Redundant Joints

Because our convenient method directly provides only the extension cylinder and joint motions for the cost functions, the method we propose includes the application of inverse kinematics to resolve the motions of the non-redundant joints. These motions are included in the vector:

$${\mathbf{x}}_k^{\left(\mathrm{r}\right)}={\left[\begin{array}{ccc}{\mathbf{q}}_k^{\left(\mathrm{r}\right)}& {\dot{\mathbf{q}}}_k^{\left(\mathrm{r}\right)}& {\ddot{\mathbf{q}}}_k^{\left(\mathrm{r}\right)}\end{array}\right]}^{\mathrm{T}}$$

(32)

For example, the positions of the base, lift and tilt joints (the non-redundant joints) are denoted with ${\mathbf{q}}_k^{\left(\mathrm{r}\right)}={\begin{array}{c}\left[q_{1,k}{q}_{2,k}{q}_{3,k}\right]\end{array}}^{\mathrm{T}}$ at time stage k. Because the extension joint’s motion is known, these non-redundant joint positions can be resolved directly by inverting Equations (4)–(6), which yields [27]:

$${\mathbf{q}}_k^{\left(\mathrm{r}\right)}={\left[\begin{array}{ccc}arctan(z_{\mathrm{w}\mathrm{d},\mathrm{k}}/x_{\mathrm{w}\mathrm{d},\mathrm{k}})& h_2& h_3\end{array}\right]}^{\mathrm{T}}$$

(33)

Here, $z_{\mathrm{w}\mathrm{d},k}$ and $x_{\mathrm{w}\mathrm{d},k}$ denote the components of the end-effector’s desired position ${\mathbf{x}}_{\mathrm{td}}$ at time stage k in the direction of the z- and x-axes, respectively, in the global frame. The lengthy inverse kinematic functions $h_2$ and $h_3$ are omitted for brevity.

At the velocity level, the approach requires that we subtract the contribution of the extension joint from the desired end-effector velocity ${\dot{\mathbf{x}}}_{\mathrm{td},k}$. Then, the joint velocities in ${\dot{\mathbf{q}}}_k^{\left(\mathrm{r}\right)}$ can be obtained by inverting Equation (7) as in:

$$\begin{array}{c}\hfill {\dot{\mathbf{q}}}_k^{\left(\mathrm{r}\right)}={\left[\begin{array}{ccc}{\mathbf{J}}_1\left({\mathbf{q}}_k\right)& {\mathbf{J}}_2\left({\mathbf{q}}_k\right)& {\mathbf{J}}_3\left({\mathbf{q}}_k\right)\end{array}\right]}^{-1}\left({\dot{\mathbf{x}}}_{\mathrm{td},k}-{\mathbf{J}}_4\left({\mathbf{q}}_k\right){\dot{q}}_{4,k}\right)\end{array}$$

(34)

because the partial Jacobian inverse matrix ${\left[{\mathbf{J}}_1\left({\mathbf{q}}_k\right){\mathbf{J}}_2\left({\mathbf{q}}_k\right){\mathbf{J}}_3\left({\mathbf{q}}_k\right)\right]}^{-1}$ is non-singular. This matrix can be symbolically computed in advance to possibly decrease the computational complexity.

Similarly, accelerations of the non-redundant joint can be resolved by inverting Equation (8):

$$\begin{array}{c}\hfill {\ddot{\mathbf{q}}}_k^{\left(\mathrm{r}\right)}={\left[\begin{array}{c}{\mathbf{J}}_1\left({\mathbf{q}}_k\right){\mathbf{J}}_2\left({\mathbf{q}}_k\right){\mathbf{J}}_3\left({\mathbf{q}}_k\right)\end{array}\right]}^{-1}\left({\ddot{\mathbf{x}}}_{\mathrm{td},k}-\dot{\mathbf{J}}({\mathbf{q}}_k,{\dot{\mathbf{q}}}_k){\dot{\mathbf{q}}}_k-{\mathbf{J}}_4\left({\mathbf{q}}_k\right){\ddot{q}}_4\right)\end{array}$$

(35)

where the right-hand term in the parentheses contains the Cartesian acceleration required from the non-redundant joints to maintain the desired end-effector acceleration ${\ddot{\mathbf{x}}}_{\mathrm{td},k}$. The Jacobian time derivative $\dot{\mathbf{J}}(\mathbf{q},\dot{\mathbf{q}})$ can be symbolically computed in advance using Equation (9). We are now in a position to apply the cylinder coordinate transformations presented to obtain the necessary motions of the cylinders of the non-redundant joints for the proposed cost functions.

The proposed method optimises the dynamical motions of all the cylinders. However, our method reduces computational complexity by only computing the necessary extension cylinder motion to optimise the dynamical motions of all the cylinders. This is possible because the motion of the extension cylinder in the 4-DOF manipulator determines the motion of the base, lift and tilt cylinders along a certain end-effector path. We could have also chosen to directly optimise the motion of the lift or tilt cylinder and use the motion of that cylinder to compute the required optimal dynamical motions of the base, lift/tilt and extension joints for the cost functions.

4.2. Proposed Cost Function: Power Consumption of Load-Sensing System

We first propose the power consumption the LS system as the principal cost in view of Equation (29) to be minimised. The cost function of the LS system, at time stage k, can be written as:

$$P_{\mathrm{h},k}=\frac{p_{\mathrm{s},k}{Q}_{\mathrm{p},k}}{{\eta }_k}$$

(36)

where $p_{\mathrm{s},k}$ denotes the pump supply pressure, $Q_{\mathrm{p},k}$ is the pump flow rate and ${\eta }_k$ is the total efficiency of the hydraulic pump and the driving motor, which depends on the operating point, i.e., the displacement, supply pressure and rotational speed of the pump. Thus, the energy consumption of the LS system can be optimised by minimising the hydraulic output power or by increasing the driving system’s efficiency. Since pump and motor efficiency data is not usually available, without loss of generality, we assume the total efficiency is constant.

The pump flow rate, which is a positive function without energy recuperation, can be approximated at time stage k by summing the cylinder flow rates as in:

$$\begin{array}{cc}\hfill Q_{\mathrm{p},k}=& \sum _{i=1}^4{v}_{i,k}\left[A_{\mathrm{A}i}H\left(v_{i,k}\right)-A_{\mathrm{B}i}H(-v_{i,k})\right]\hfill \end{array}$$

(37)

where $v_{i,k}$ is the velocity of cylinder i and $H\left(v_{i,k}\right)$ denotes the piecewise Heaviside step function. The piston-side area is denoted with $A_{\mathrm{A}i}$ and the rod-side area with $A_{\mathrm{B}i}$. For simplicity, the hydraulic fluid is assumed to be incompressible. The rod-side area is negatively signed to enforce a positive pump flow rate when the cylinder is retracting and has negative velocity. The Heaviside step function can be defined using:

$$H\left(v_{i,k}\right)=\left\{\begin{array}{cc}0\hfill & \mathrm{if}{v}_{i,k}<0\hfill \\ \frac{1}{2}\hfill & \mathrm{if}{v}_{i,k}=0\hfill \\ 1\hfill & \mathrm{if}{v}_{i,k}>0\hfill \end{array}\right.$$

(38)

which introduces a discontinuity in the objective function. Because having a differentiable objective function is desirable, the discontinuous Heaviside step function can be approximated by a differentiable expression, after some mathematical manipulations:

$$H\left(v_{i,k}\right)\approx \widehat{H}\left(v_{i,k}\right)=\frac{1}{1+e^{-2s{v}_{i,k}}}$$

(39)

where s is a real positive number, such as 20,000.

The pump supply pressure, according to the LS system’s operating principle, can be calculated at steady-state at time stage k by using the maximum driven actuator pressure:

$$p_{\mathrm{s},k}=max\left\{p_{1,k},p_{2,k},p_{3,k},p_{4,k}\right\}+\Delta p_{\mathrm{LS}}$$

(40)

where $p_{1,k}$, $p_{2,k}$, $p_{3,k}$ and $p_{4,k}$ denote the chamber pressures of the base, lift, tilt and extension cylinders, respectively. The LS pressure margin $\Delta p_{\mathrm{LS}}$ is roughly 3 MPa in the test-bed.

The chamber pressure $p_{i,k}$ of cylinder i is taken from the rod- or piston-side depending on the cylinder’s direction of motion. The pressure is zero if the actuator is not driven. Hence, the chamber pressure $p_{i,k}$ at time stage k can be estimated by using:

$$p_i=\left(p_{\mathrm{ssA},i}H(v_i-\epsilon )+p_{\mathrm{ssB},i}H(-v_i-\epsilon )\right)H\left(F_i{v}_i\right)$$

(41)

where the use of the minor positive constant $\epsilon$ in the Heaviside step functions assures zero chamber pressure when the cylinder is not driven and its velocity $v_i\in [-\epsilon ,\epsilon ]$. Time indices were omitted for clarity. In addition, multiplying by the Heaviside value of the actuator power ensures a non-negative actuator pressure when the cylinder load is over-running. In this over-running load case, the cylinder velocity $v_i$ and the cylinder pressure force $F_i$ have opposite signs. The steady-state pressure equations are chosen according to the spool type used either from Equations (23)–(24) or Equations (27)–(28).

A discontinuous maximum function has been used in Equation (40). This function can be approximated by applying the equation:

$$max(p_1,p_2)=0.5\left(p_1+p_2+|p_1-p_2|\right)$$

(42)

thrice; i.e., $max\{p_{1,k},p_{2,k},p_{3,k},p_{4,k}\}$ can be computed as $max(max(max(p_{1,k},p_{2,k}),p_{3,k}),p_{4,k})$. The absolute function, in turn, can be approximated as the square root of ${(p_1-p_2)}^2+ϵ$. This $ϵ$ should be a minor positive constant, such as ${10}^{-6}$, to improve the accuracy of the approximation. As a result of these approximations, the LS cost is desirably continuously differentiable, which allows us to more easily employ a standard optimal control solver for performance comparison. The pump flow rate and supply pressure modelling results achieved with the test-bed can be found from Appendix C.

4.3. Proposed Cost Function: Power Consumption of Constant-Pressure System

In our second application, we minimise the energy consumption of the CP system. The cost function of the CP system can be written as:

$$P_{\mathrm{h},k}=\frac{p_{\mathrm{s}}{Q}_{\mathrm{p},k}}{{\eta }_k}$$

(43)

where $p_{\mathrm{s}}$ denotes the constant supply pressure and the pump flow rate $Q_{\mathrm{p},k}$ is given by Equation (37). Because of the constant pressure, the energy consumption of the CP system can generally be optimised load-independently without knowledge of manipulator dynamics through minimisation of the pump flow rate. By minimising the pump flow rate, the total volumes displaced by the hydraulic cylinders in the solution end up being minimised. Because the piston and piston-rod-side areas obviously affect the cylinder volumes displaced, minimising the energy consumption of the CP system does not equal to the actuator velocity norm minimisation, e.g., in [10]. The variation of cylinder areas among the manipulator’s hydraulic cylinders, including the different piston and rod-side areas of these conventionally asymmetric cylinders, render the pump flow rate minimisation a highly tractable objective.

4.4. Actuator Limits

Because of the inherent physical limitations of the manipulator’s cylinders, the cylinders’ motions are limited. For this reason, we have introduced a mechanism through the optimal control problem’s term ${\ell }_k({\mathbf{q}}_k,{\dot{\mathbf{q}}}_k,{\ddot{\mathbf{q}}}_k)$ to impose actuator constraints at the position, velocity and acceleration levels by issuing a significant penalty compared to the magnitude of the principal cost if the limits are violated.

Since the position limits depend on the cylinder stroke, they can be expressed effectively either in the joint or actuator space. In contrast, however, the velocity limits should be expressed in the actuator space because these limits cannot be expressed optimally from the perspective of the hydraulic system in the joint space. Based on the influential dimensions of the hydraulic pump and cylinders and the cylinders’ non-constant torque arms these velocity limits are constant only in the actuator space to ensure that the total fluid flow rates for the actuator do not exceed the pump’s maximum flow output. Assuming fluid incompressibility in these computations is the most convenient approach. Notice that by assuming constant velocity limits, the joint velocity limits vary if the torque arm is non-constant (see Equation (11)). Only if a particular cylinder’s torque arm is constant can the corresponding velocity limit be expressed optimally in the joint space. Hence, in the test-bed case, only the velocity limits for base and extension cylinders can be expressed in either space.

Lastly, we choose constant acceleration limits for the cylinders for simplicity. Consequently, the joint acceleration limits vary (see Equation (12)). Selecting the proper acceleration limit magnitudes are important from at least two perspectives: firstly, to limit the occurrence of higher accelerations that might lower the manipulator wear life and secondly, to assure higher controller tracking performance. The actuator limits used in the experiments can be found in

Table 3. Actuator/joint limits of the test-bed manipulator.

i	Joint Pos.	Cyl. Vel. [m/s]	Cyl. Acc. [m${}^2$/s]
1	$[-\pi ,\pi ]$ rad	$[-0.10,0.10]$	$[-0.4,0.4]$
2	$[-0.425,1.097]$ rad	$[-0.15,0.15]$	$[-0.4,0.4]$
3	$[-2.75,-0.25]$ rad	$[-0.15,0.15]$	$[-0.4,0.4]$
4	$[0,2\times 1.85]$ m	$[-0.40,0.40]$	$[-0.4,0.4]$

.

5. Dynamic Programming Solution

The complex optimal control problem of the hydraulic manipulator is solved to global optimality using the discrete-time DP algorithm (see Algorithm 1 [31,32,33]). The main advantage of the algorithm is that it can produce this solution considerably more effectively than a brute-force method, whose computational cost grows exponentially with the time stages N. However, the main disadvantage of the algorithm is that its computational cost grows exponentially with the states and controls due to the curse of dimensionality. Hence, our low-dimensional optimal control problem that we have reduced to the bare essentials makes this global DP algorithm applicable. Recall that we have reduced our optimal control problem size to only two states and one control.

The effectiveness of the DP algorithm over the brute-force method results from the principle of optimality, which states that “an optimal policy has the property that whatever the initial state and initial decision (control) are, the remaining decisions must constitute an optimal policy with regard to the state resulting from the first decision” [31]. The Bellman equation, which characterises this principle of optimality, can be written at time stage k, in view of our objective Equation (29), as:

$$\begin{array}{c}\hfill J_k^{☆}\left({\overline{\mathbf{x}}}_k\right)=\underset{{\overline{u}}_k\in \mathbb{U}}{min}\left\{P_{\mathrm{h},k}\left({\mathbf{q}}_k,{\dot{\mathbf{q}}}_k,{\ddot{\mathbf{q}}}_k\right)+{\ell }_k\left({\mathbf{q}}_k,{\dot{\mathbf{q}}}_k,{\ddot{\mathbf{q}}}_k\right)+J_{k+1}^{☆}\left({\overline{\mathbf{F}}}_k({\overline{\mathbf{x}}}_k,{\overline{u}}_k)\right)\right\}\end{array}$$

(44)

where the sum $P_{\mathrm{h},k}+{\ell }_k$ represents the running cost and the vector ${\overline{\mathbf{F}}}_k$ represents the system dynamics (see Equations (30)–(31)). The minimum cost-to-go $J_k^{☆}$ is the minimum sum of the running costs from stage $k\in \{0,1,\dots ,N-1\}$ to the final stage at $k=N-1$. This optimal cost-to-go $J_k^{☆}\left({\overline{\mathbf{x}}}_k\right)$ is defined for each state vector combination such that ${\overline{\mathbf{x}}}_k\in \mathbf{𝕏}=\{{\mathbb{X}}_1\times {\mathbb{X}}_2\}$ and minimised over all possible controls at stage k such that ${\overline{u}}_k\in \mathbb{U}$. Here, ${\mathbb{X}}_i=\{x_i^{\left(1\right)},x_i^{\left(2\right)},\dots ,x_i^{\left(N_{x_i}\right)}\}$ and $\mathbb{U}=\{u^{\left(1\right)},u^{\left(2\right)},\dots ,u^{\left(N_u\right)}\}$ denote the discrete-state and control sets, respectively. The variables $N_{{\mathrm{x}}_i}$ and $N_{\mathrm{u}}$, in turn, determine the density of the state and control grids, respectively. Sufficiently dense grids are defined to reach the global optimum. These grid sizes, and suitable time-discretisation, can be found through trial and error, taking note that the computational and memory requirements of the algorithm scale linearly with the density of the grids.

The algorithm’s implementation shown in Algorithm 1 entails two consecutive phases: a backwards phase followed by a forward simulation. During the computationally heavier backwards phase, the optimal controls are resolved at each position in the state grid while the algorithm iterates backwards in time. In greater detail, at any given time stage k, the system is simulated (line 5), subsequent computations are performed to evaluate the cost function (lines 6–8), the minimum cost-to-go $J_k^{☆}\left({\overline{\mathbf{x}}}_{k+1}\right)$ computed at the previous iteration is interpolated because the system dynamics have likely transformed into a state unspecified in the grid (line 9) and the cost-to-go $J_k\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{u}}}_k\right)$ is evaluated for all the states and controls by utilising the minimum cost-to-go $J_{k+1}^{☆}\left({\overline{\mathbf{x}}}_{k+1}\right)$ (line 10). From the cost-to-go evaluated at each control, the minimum cost-to-go $J_k^{☆}\left({\overline{\mathbf{x}}}_k\right)$ is then stored for the next iteration (line 12) together with the optimal control $u_k^{☆}\left({\overline{\mathbf{x}}}_k\right)$, which minimised the cost-to-go (line 13). The backwards phase iterates in this manner until it reaches the initial time stage. At this point, with our system of two states and one control, the algorithm has produced a two-dimensional state grid at each time stage, and the optimal controls have been stored at each position in the state grid. The algorithm then simulates in the forward direction from the initial time stage and state by following the grid of the optimal controls produced in the backwards phase. Interpolation of the optimal control is used during the forward simulation because is it is unlikely that the optimal controls were evaluated at the optimal states. When the simulation has been completed, the algorithm has produced the optimal control sequence and state trajectories as outputs. In addition to producing the global solution, this algorithm has the significant advantage that different initial configurations of the manipulator can be considered with just a forward simulation. This advantage is in contrast to most other optimisation algorithms where new optimisations from scratch are required.

Algorithm 1 Discrete-time dynamic programming algorithm

1. Backward phase with interpolation Require: Sets $\mathbf{𝕏}$ and $\mathbb{U}$ 1: Initialisation: set final cost-to-go $J_N^{☆}\left({\overline{\mathbf{x}}}_N\right)$ for all ${\overline{\mathbf{x}}}_N\in \mathbf{𝕏}$ 2: for $k=N-1$ to 0 do 3:   for all ${\overline{\mathbf{x}}}_k\in \mathbf{𝕏}$ do 4:     for all ${\overline{u}}_k\in \mathbb{U}$ do 5:      ${\overline{\mathbf{x}}}_{k+1}={\overline{\mathbf{F}}}_k({\overline{\mathbf{x}}}_k,{\overline{u}}_k)$, see Equations (30) and (31). 6:      Determine $q_{4,k},{\dot{q}}_{4,k}\hspace{0.17em}\mathrm{and}\hspace{0.17em}{\ddot{q}}_{4,k}$ using Equations (10)–(12). 7:      Complete ${\mathbf{q}}_k,{\dot{\mathbf{q}}}_{k}and{\ddot{\mathbf{q}}}_k$ using Equations (33)–(35). 8:      Compute base, lift and tilt cyl. coordinates using Equations (10)–(12). 9:      Interpolate $J_{k+1}^{☆}\left({\overline{\mathbf{x}}}_{k+1}\right)$ 10:      $J_k\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{u}}}_k\right)=P_{\mathrm{h},k}({\mathbf{q}}_k,{\dot{\mathbf{q}}}_k,{\ddot{\mathbf{q}}}_k)+{\ell }_k({\mathbf{q}}_k,{\dot{\mathbf{q}}}_k,{\ddot{\mathbf{q}}}_k)+J_{k+1}^{☆}\left({\overline{\mathbf{x}}}_{k+1}\right)$ 11:    end for 12:    $J_k^{☆}\left({\overline{\mathbf{x}}}_k\right)={min}_{{\overline{u}}_k\in \mathbb{U}}{J}_k\left({\overline{\mathbf{x}}}_k,{\overline{u}}_k\right)$ 13:    $u_k^{☆}\left({\overline{\mathbf{x}}}_k\right)={\mathrm{argmin}}_{{\overline{u}}_k\in \mathbb{U}}{J}_k\left({\overline{\mathbf{x}}}_k,{\overline{u}}_k\right)$ 14:   end for 15: end for 2. Forward simulation with interpolation Require: $u_k^{☆}\left({\overline{\mathbf{x}}}_k\right)$ for all $k\in \{0,1,\dots ,N-1\}$ and ${\overline{\mathbf{x}}}_k\in \mathbf{𝕏}$ 1: Initialisation: set ${\overline{\mathbf{x}}}_0$ 2: for $k=0$ to $N-1$ do 3:   Interpolate $u_k^{☆}\left({\overline{\mathbf{x}}}_k\right)$ 4:   ${\overline{\mathbf{x}}}_{k+1}={\overline{\mathbf{F}}}_k\left({\overline{\mathbf{x}}}_k,u_k^{☆}\left({\overline{\mathbf{x}}}_k\right)\right)$ 5: end for 6: Determine $q_{4,k},{\dot{q}}_{4,k}\mathrm{and}{\ddot{q}}_{4,k}$ using Equations (10)–(12). 7: Complete ${\mathbf{q}}_k,{\dot{\mathbf{q}}}_k\mathrm{and}{\ddot{\mathbf{q}}}_k$ using Equations (33)–(35).

A demonstratively simple example of the DP algorithm has been presented, e.g., by Nurmi and Mattila [21]. The utilised general Matlab-coded DP algorithm can be found from Sundstrom and Guzzella [32]. Different DP algorithm varieties have been compared by van Berkel et al. [34] and Böhme and Frank [33], and a convergence proof of discrete-time DP has been presented by Shin and McKay [35].

6. Experimental Results

Before diving head-first into the experimental results from the real hydraulic manipulator introduced in Section 2, we describe the testing methodology. Firstly, the optimal control solutions compared consist of:

(1)

the proposed LS cost (36), which is minimised by using the DP algorithm and named DP (LS) according to the naming convention “method (minimised cost)”;

(2)

the proposed LS cost, which is minimised by using the DIDO algorithm [22] and named DIDO (LS);

(3)

the proposed CP cost (43), which is minimised according to DP (CP);

(4)

the proposed CP cost, which is minimised according to DIDO (CP);

(5)

a positive actuator energy cost, which is solved according to DP (A. En.) and which minimises $P_{\mathrm{h},k}={\sum }_{i=1}^4{F}_{i,k}{v}_{i,k}$ for positive work, $F_{i,k}{v}_{i,k}>0$, and

(6)

an extension cylinder use cost, which is solved according to DP (ext.) and which minimises $P_{\mathrm{h},k}=v_{4,k}^2$.

Secondly, the instantaneously optimal solutions included in the comparison are

(7)

the joint-limited pseudo-inverse solution by Flacco et al. [20] named P-inv. (j.), which minimises ${\dot{\mathbf{q}}}^{\mathrm{T}}\dot{\mathbf{q}}$, and

(8)

an actuator-limited pseudo-inverse solution named P-inv. (a.), which was derived by combining the solutions of Beiner and Mattila [10] and Flacco et al. [20] to minimise ${\mathbf{v}}^{\mathrm{T}}\mathbf{v}$.

The DP solution (A. En.) is formulated with positive actuator work only to obtain the best attainable performance in the test-bed, which did not have the possibility of recovering energy. The solution for minimising use of the extension cylinder is included in the comparison to evaluate whether the solution lowers energy consumption, presumably, by reducing the lift and tilt cylinder forces required.

The DP algorithm was parametrised by dividing the state and control grid into 356 × 101 × 151 ($N_{{\mathrm{x}}_1}$ × $N_{{\mathrm{x}}_2}$ × $N_{\mathrm{u}}$) discrete values. With this grid size and by using a time discretisation of 0.05 s, we could nearly max out the 8 GB of RAM available. In the evaluation, the performance of the DP algorithm is compared to a “bootstrapped” DIDO solver by employing the proposed problem formulation. This DIDO algorithm is a commercial, pseudospectral and Matlab-based optimal control solver that has been extensively used in the solution for aerospace optimal control problems.

The energy consumption of the LS and CP systems was computed as the performance indices to be compared by measuring the pump supply pressure and the flow rate along three kinds of end-effector paths. The end-effector was controlled without payload along these paths. Firstly, we present the energy consumption in the horizontal and vertical point-to-point paths where the end-effector moves from point A to point B without a payload in a straight line (see Figure 9a). Jazar’s quintic rest-to-rest polynomial algorithm [36] was used to generate the Cartesian motion profile between these points. The horizontal path was driven in 10 s, requiring a maximum end-effector velocity of 0.60 m/s along the $x_{\mathrm{w}}$-axis. The vertical path was driven in 6 s, and it required a higher maximum end-effector velocity of 0.95 m/s along the $y_{\mathrm{w}}$-axis. Secondly, we consider half of the log-loading cycle where the end-effector starts from the load space (or trailer) of the forwarder at point A and moves the empty log grapple back to the logs on the ground at point B (see Figure 9b). A cubic spline was used to generate the desired end-effector path via points in this case (see e.g., Biagiotti and Melchiorri [37]). This path was driven in 15 s, and the maximum end-effector velocity was 0.75 m/s in the $x_{\mathrm{w}}{y}_{\mathrm{w}}{z}_{\mathrm{w}}$-space.

The paths were driven three times using each redundancy resolution from (1) through (8), and the results were averaged. Two initial configurations of the manipulator at the beginning of each path were also considered: firstly, based on the smallest feasible extension cylinder length and subsequently based on the middle point between the smallest and largest feasible extension cylinder lengths. The supply pressure of the CP system was set to 20 MPa in the experiments, and the motion controller used is discussed in Appendix A.

6.1. Horizontal Path

The experimentally obtained energy consumption of the LS system of each redundancy resolution has been sorted by magnitude in Figure 10a. These results demonstrate that the pseudo-inverses clearly have a sub-optimal performance by yielding 25% greater energy consumption than the best proposed hydraulic-energy optimal solution. The DIDO and DP solutions produce comparable results. In this path, minimising the energy consumption of the actuators is also similar to minimising the energy consumption of the LS system. Another important observation is the adequate performance of the CP solution, considering that it is based on a simple model. In addition, because the extension joint use minimising cost yields a poor solution, minimising the use of the extension joint appears to be impractical for our goal of minimising the hydraulic energy consumption.

As can be seen in Figure 10a, the CP system results are comparable to the LS results, but the main difference is the greater energy consumption. The pseudo-inverses, again, have an inferior performance from the perspective of the energy consumption of the CP system, and they require approximately 30% greater energy use than do the hydraulic energy-optimal methods in this case. Inferiority of the pseudo-inverses reinforces the notion that minimising the velocities of the actuators does not minimise the energy consumption of the CP system. In addition, the best solutions here consume the same amount of CP energy as do the pseudo-inverses in the LS case with the same path.

Two of the redundancy resolutions’ joint trajectories and controller performances have been illustrated in Figure 10b,c. The DP joint references, as shown there, were always low-pass filtered (cutoff frequency 5 Hz) to smooth out the trajectories’ slightly jagged edges. This filtering caused no apparent error in the desired end-effector path. The results with the second initial configuration can be found in Appendix D (see Figure A4).

6.2. Vertical Path

Continuing with the results of the vertical path, we firstly compare the sorted LS energy-consumption presented in Figure 11a. As previously, the energy-optimised solutions for the hydraulic manipulators are superior compared to the pseudo-inverses, which require almost twice the energy consumption. However, the DP energy-optimised solutions for hydraulic manipulators have not reached a global optimum because the DIDO method clearly fares better in this case. The DIDO solutions themselves are practically the same. This path also demonstrates that minimising the energy consumption of the actuators does not always minimise the energy consumption of the LS system. The proposed solutions minimising the CP cost, again, have an adequate performance.

The energy-consumptions rates of the CP system presented in Figure 11a are similar, but the energy-consumption rates are approximately double compared to the LS case. The pseudo-inverses, most evidently, have an inferior performance compared to the proposed solutions, which, in contrast, require the lowest energy consumption.

A few optimal and sub-optimal joint trajectories and controller performances are illustrated in Figure 11b,c. Compared to the previous path, the controller tracking performances are slightly worse because the end-effector moves faster. The energy consumption with the second initial configuration is presented in Figure A5 in Appendix D. These energy consumption rates demonstrate that the initial manipulator configuration has a significant impact on the energy consumption. Figure A5 also shows that the algorithms, or DP (CP) specifically in this case, can sometimes yield an average result. We suspect that this result is related to a combination of the DP algorithm’s grid size, closed-loop dynamics and modelling errors.

6.3. Half-Cycle

Finally, we have the “half-cycle” energy consumption. The energy-consumption rates of the LS system presented in Figure 12a demonstrate that the redundancy resolutions’ perform fairly similarly in this case. The minimisation of the actuator energy is also optimal in this case. The joint-limited pseudo-inverse, surprisingly, produces an energy-consumption rate comparable to the consumption obtained using hydraulic energy-optimised redundancy resolutions. Because the simulations (not presented in this paper) predicted a greater difference between these solutions, the smaller difference originates from modelling errors and the behaviour of the closed-loop control. The energy-consumption rates for the CP system presented in Figure 12a differ from the LS results in that the pseudo-inverses are inferior.

Selected optimal and sub-optimal joint trajectories and controller performances from the half-cycle are illustrated in Figure 12b,c. Furthermore, the controller performances of each test case were evaluated according to the ISO 9283 standard, which contains a standardised general method for computing the maximum end-effector path tracking error (see

Table 4. Maximum end-effector path tracking errors in each redundancy resolution computed according to ISO 9283 standard [38].

Path	Solution	Track. Err. (m)
		Load-Sensing System	Constant-Pressure System
Horizontal	DP (CP)	0.088	0.070
	DP (LS)	0.076	0.073
	DP (A. En.)	0.064	0.064
	DP (ext.)	0.146	0.140
	P-inv. (a.)	0.127	0.069
	P-inv. (j.)	0.057	0.052
	DIDO (CP)	0.064	0.041
	DIDO (LS)	0.063	0.049
Vertical	DP (CP)	0.098	0.091
	DP (LS)	0.093	0.087
	DP (A. En.)	0.076	0.085
	DP (ext.)	0.085	0.086
	P-inv. (a.)	0.108	0.095
	P-inv. (j.)	0.098	0.101
	DIDO (CP)	0.077	0.080
	DIDO (LS)	0.079	0.081
Half-cycle	DP (CP)	0.346	0.129
	DP (LS)	0.265	0.117
	DP (A. En.)	0.232	0.111
	DP (ext.)	0.210	0.155
	P-inv. (a.)	0.200	0.219
	P-inv. (j.)	0.343	0.207
	DIDO (CP)	0.277	0.148
	DIDO (LS)	0.245	0.179

). The controller performances are superior when the CP hydraulic system is used, particularly over the half-cycle path. This is not surprising because the LS system’s actuator flow dynamics are slower due to the non-constant pressure level. The tracking errors, in general, are sufficiently insignificant and consistent for evaluating the proposed solution. Again, the result with the second initial configuration can be found in Appendix D (Figure A6).

As a conclusion of the evaluation, the pseudo-inverses, in particular, and actuator energy solutions, in some cases, were inferior compared to the proposed hydraulic energy-optimised redundancy resolutions. These inferior solutions should be used with caution when minimising hydraulic energy consumption. The DIDO solutions were, for the most part, comparable to the DP solutions. The CP optimal solutions, in turn, were often comparable to the LS optimal solutions. The vertical and horizontal paths demonstrated the greatest energy savings, but energy savings were also attainable in the half-cycle when energy-optimised redundancy resolution was used. A solution containing both halves of the log-loading cycle could be obtained from a combination of a minimum-time solution for the log-loading part and the energy-optimal solution proposed for this return part. The results were obtained with the assumption of constant pumping efficiency. In the conditions of the experiments, in which the pump flow rate was typically less than 25% of the maximum flow rate of the pump, it might have been advantageous to incorporate the variations in pumping efficiency in the optimisation to improve the results. However, the satisfactory results of our paper suggest that it might not be worthwhile to consider the efficiency map of the pump in the optimisation.

7. Conclusions and Future Work

This paper addressed the open problem of global energy-optimal redundancy resolution of 4-DOF hydraulic manipulators in robotic control applications, in which the highly repetitive manipulator movements are automated to improve work productivity and operator workload circumstances. To ensure the hydraulic energy-optimality of the redundancy resolution, the solution was formulated at the hydraulic system level by proposing effective cost functions for the commonly used LS and CP hydraulic systems. Furthermore, actuator constraints at the position, velocity and acceleration levels were imposed in the solution. To reduce the computational complexity and problem dimensionality, firstly, only the motion of the redundant extension cylinder was directly optimised. Secondly, hydraulic system properties were incorporated through steady-state equations in the proposed cost functions. The non-redundant actuator motions were resolved via non-singular inverse kinematics to obtain the full manipulator motion state for the cost functions. The global DP algorithm was applicable to our optimal control problem formulation due to the conveniently reduced problem dimensionality. A pseudocode implementation of the DP algorithm used in the optimal control problem solution was presented to discuss the properties of the algorithm and to showcase its relatively simple implementation. Field experiments that were performed on a 4-DOF forestry manipulator demonstrated the significant advantage of the new global optimal control problem solution. Around 15–30% greater hydraulic energy consumption was observed with the conventional solutions in the LS and CP systems. These experiments also showed that the CP system’s cost function for basic energy consumption, most often, also led to minimised energy consumption by the LS system. Overall, these experimental results demonstrate energy-optimal redundancy resolution at the hydraulic system level in the prospective automated robotic applications. The existing practices in redundancy resolution of hydraulic manipulators are, therefore, significantly expanded. Experimental evaluation of the solution performance with an end-effector load was left for future work to limit the study’s scope.