Sub-optimal Control of Autonomous Wheel loader with Approximate Dynamic Programming

TL;DR

This paper presents a method for controlling an autonomous wheel loader during short loading cycles by optimizing mode switching times using approximate dynamic programming, based on a detailed vehicle and engine model.

Contribution

It introduces a control approach for wheel loaders that optimizes switching times between modes using approximate dynamic programming, considering a fixed mode sequence and known optimal path.

Findings

Simulation results demonstrate the effectiveness of the control strategy.

The approach successfully optimizes switching times for improved performance.

Abstract

Optimal control of wheel loaders in short loading cycles is studied in this paper. For modeling the wheel loader, the data from a validated diesel engine model is used to find a control oriented mean value engine model. The driveline is modeled as a switched system with three constant gear ratios (modes) of $-60$ for backwarding, $60$ for forwarding, and zero for stopping. With these three modes, the sequence of active modes in a short loading cycle is fixed as backwarding, stopping, forwarding, and stopping. For the control part, it is assumed that the optimal path is known a priori. Given the mode sequence, the control objective is finding the optimal switching time instants between the modes while the wheel loader tracks the optimal path. To solve the optimal control problem, approximate dynamic programming is used. Simulation results are provided to show the effectiveness of the solution.

Tables

Table 1. TABLE I: Nomenclature

${\omega }_e$	rotational speed of engine	$\alpha$	angle between lift force and lift arm
$J_e$	engine moment of inertia	$\theta$	lift arm angle
$T_e$	engine torque	$T_{buc}$	torque due to the bucket load
$p_{im}$	intake manifold pressure	$T_{arm,w}$	torque due to boom weight
$u_f$	mass of fuel per engine cycle	$I_{boom}$	boom moment of inertia
$\omega$	angular velocity of the boom	$F_w$	wheels’ traction force
$V$	vehicle velocity	$F_{roll}$	rolling resistance force
$u_p$	pressure in the lift cylinder	$M_{tot}$	mass of WL and load
$u_s$	time derivative of steering angle	$X$	position of WL
${\tau }_p$	time constant in $p_{im}$	$Y$	position of WL
$p_{stat}$	stationary intake manifold pressure	$\beta$	heading angle
$F_{cyl}$	lifting force	$R_{turn}$	turning radius
$R$	a parameter in the boom geometry	$\delta$	steering angle

Table 2. TABLE II: Parameters of WL model

Parameter	Description
$J_e$ = 0.43	engine inertia
$r_1$ = 2	lift arm dimensions
$r_2$ = 2	lift arm dimensions
${\theta }_1$ = pi/6	bend angle of the lift arm
$n_{lc}$ = 2	number of lift cylinders
$A_{lc}$ = 0.0284	lift cylinder cross section area
$y_g$ = 2.13	height between body and lift arm hinge
${\eta }_{lift}$ = 0.5	lift system efficiency
$y_{off}$ = 0.5	boom hinge offset
$M_{tot}$ = 31330	mass of WL and load
${\mu }_{roll}$ = 0.03	rolling resistant factor
$F_{roll}$ = $9.81{\mu }_{roll}{M}_{tot}$	rolling resistance
$L$ = 3.7	WL wheel base
${\rho }_f$ = 832	fuel density
$R_w$ = 0.3175	wheel radius
$T_1$ = 5	constant for $P_{pump}$
$T_2$ = -2.5	constant for $P_{pump}$
${\eta }_{gb}$ = 0.9	gearbox efficiency
$C_{st}$ = ${10}^5$	constant for steering power
$F_{buc}$ = 0.0981 $M_{tot}$	bucket and load weight
$F_{arm}$ = $F_{buc}$	arm weight
$I_{boom}$ = 1200	moment of inertia for the boom

Equations

$${\dot{\omega }}_e$$

$${\dot{p}}_{im}$$

$$\dot{\omega }$$

$$\dot{\theta }$$

$$\dot{V}$$

$$\dot{X}$$

$$\dot{Y}$$

$$\dot{\beta }$$

$$\dot{\delta }$$

$$\begin{array}{r}{\hat{T}}_e({\omega }_e,p_{im},u_f)=W_{T_e}^T[1,{\omega }_e,p_{im},u_f{]}^T\end{array}$$

$$\begin{array}{r}{P}_{e,load}=P_{lift}(\omega ,u_p)+P_{steering}(u_s)+P_{traction}({\omega }_e,V)\end{array}$$

$$P_{lift}(\omega,u_{p})=\max\Big{(}0,\frac{Q(\omega)u_{p}}{\eta_{lift}}\Big{)}$$

$$\begin{array}{r}{P}_{steering}=C_{st}{u}_s^2\end{array}$$

$$\begin{array}{r}{P}_{traction}(\omega ,V)=P_{pump}({\omega }_e,V,\gamma )|sign(\gamma )|\end{array}$$

$$\begin{array}{r}{\tau }_p({\omega }_e)={\tau }_1{\omega }_e+{\tau }_2\end{array}$$

$$\begin{array}{r}{p}_{stat}=a_1{\omega }_e+a_2{u}_f+a_3\end{array}$$

$$\begin{split}\dot{x}(t)&=\bar{f}_{v}\big{(}x(t)\big{)}+\bar{g}_{v}\big{(}x(t)\big{)}u(t),\>v\in=\{1,2,3\},\>x(0)=x_{0}\end{split}$$

$$\begin{array}{r}x(t)=[x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9{]}^T\end{array}$$

$${\dot{x}}_{10}$$

$${\dot{x}}_{11}$$

$$\begin{split}J&(x_{0})=\big{(}x(t_{f})-r(t_{f})\big{)}^{T}S\big{(}x(t_{f})-r(t_{f})\big{)}+\int_{t_{0}}^{t_{f}}\frac{1}{2}\Big{(}\big{(}x(t)-r(t)\big{)}^{T}\bar{Q}\big{(}x(t)-r(t)\big{)}+u(t)^{T}\bar{R}u(t)\Big{)}dt\end{split}$$

$$t=\begin{cases*}t_{0}+(t_{1}-t_{0})\hat{t}&if $0\leq\hat{t}<1$\\ t_{1}+(t_{f}-t_{1})(\hat{t}-1)&if $1\leq\hat{t}\leq 2$\end{cases*}$$

$$x^{\prime}(\hat{t})=\begin{cases*}\Big{(}\bar{f}_{1}\big{(}x(\hat{t})\big{)}+\bar{g}_{1}\big{(}x(\hat{t})\big{)}u(\hat{t})\Big{)}(t_{1}-t_{0})&if $0\leq\hat{t}<1$\\ \Big{(}\bar{f}_{2}\big{(}x(\hat{t})\big{)}+\bar{g}_{2}\big{(}x(\hat{t})\big{)}u(\hat{t})\Big{)}(t_{f}-t_{1})&if $1\leq\hat{t}\leq 2$\end{cases*}$$

$$\begin{split}J(x_{0})&\equiv J(t_{1},x_{0})=\big{(}x(2)-r(2)\big{)}^{T}S\big{(}x(2)-r(2)\big{)}+\int_{0}^{1}\frac{1}{2}\Big{(}\big{(}x(\hat{t})-r(\hat{t})\big{)}^{T}\bar{Q}(t_{1}-t_{0})\big{(}x(\hat{t})-r(\hat{t})\big{)}+u(\hat{t})^{T}\bar{R}(t_{1}-t_{0})u(\hat{t})\Big{)}d\hat{t}\\ &+\int_{1}^{2}\frac{1}{2}\Big{(}\big{(}x(\hat{t})-r(\hat{t})\big{)}^{T}\bar{Q}(t_{f}-t_{1})\big{(}x(\hat{t})-r(\hat{t})\big{)}+u(\hat{t})^{T}\bar{R}(t_{f}-t_{1})u(\hat{t})\Big{)}d\hat{t}\end{split}$$

$$x_{\hat{k}+1}=\begin{cases*}f_{1}\big{(}x_{\hat{k}}\big{)}+g_{1}\big{(}x_{\hat{k}}\big{)}u_{\hat{k}}&if $0\leq\hat{k}<\frac{N^{\prime}}{2}$\\ f_{2}\big{(}x_{\hat{k}}\big{)}+g_{2}\big{(}x_{\hat{k}}\big{)}u_{\hat{k}}&if $\frac{N^{\prime}}{2}\leq\hat{k}\leq N^{\prime}$\end{cases*}$$

$$\begin{split}J(t_{1},x_{0})&=\big{(}x_{N^{\prime}}-r_{N^{\prime}}\big{)}^{T}S\big{(}x_{N^{\prime}}-r_{N^{\prime}}\big{)}\\ &+\frac{1}{2}\sum_{\hat{k}=1}^{N^{\prime}/2}\Big{(}(x_{\hat{k}}-r_{\hat{k}})^{T}\bar{Q}(t_{1}-t_{0})\delta\hat{t}(x_{\hat{k}}-r_{\hat{k}})+u_{\hat{k}}^{T}\bar{R}(t_{1}-t_{0})\delta\hat{t}u_{\hat{k}}\Big{)}\\ &+\frac{1}{2}\sum_{\hat{k}=N^{\prime}/2+1}^{N^{\prime}-1}\Big{(}(x_{\hat{k}}-r_{\hat{k}})^{T}\bar{Q}(t_{f}-t_{1})\delta\hat{t}(x_{\hat{k}}-r_{\hat{k}})+u_{\hat{k}}^{T}\bar{R}(t_{f}-t_{1})\delta\hat{t}u_{\hat{k}}\Big{)}\end{split}$$

$$\begin{split}V(t_{1},x_{\hat{k}})=&\min_{u(.)}\Big{(}\frac{1}{2}(x_{\hat{k}}-r_{\hat{k}})^{T}\bar{Q}(t_{1}-t_{0})\delta\hat{t}(x_{\hat{k}}-r_{\hat{k}})+\frac{1}{2}u_{\hat{k}}^{T}\bar{R}(t_{1}-t_{0})\delta\hat{t}u_{\hat{k}}+V(t_{1},x_{\hat{k}+1})\Big{)}\end{split}$$

$$\begin{split}u^{*}_{{\hat{k}}}=\operatorname*{arg\,min}_{u(.)}\Big{(}&\frac{1}{2}(x_{\hat{k}}-r_{\hat{k}})^{T}\bar{Q}(t_{1}-t_{0})\delta\hat{t}(x_{\hat{k}}-r_{\hat{k}})+\frac{1}{2}u_{\hat{k}}^{T}\bar{R}(t_{1}-t_{0})\delta\hat{t}u_{\hat{k}}+V(t_{1},x_{\hat{k}+1})\Big{)}\end{split}$$

$$\lambda_{\hat{k}}(t_{1},x_{\hat{k}})=\begin{cases*}{\mathscr{Q}}_{1}+\frac{\partial x^{*}_{\hat{k}+1}}{\partial x_{\hat{k}}}\lambda_{\hat{k}+1}(t_{1},x^{*}_{\hat{k}+1})&if $0\leq\hat{k}<\frac{N^{\prime}}{2}$\\ {\mathscr{Q}}_{2}+\frac{\partial x^{*}_{\hat{k}+1}}{\partial x_{\hat{k}}}\lambda_{\hat{k}+1}(t_{1},x^{*}_{\hat{k}+1})&if $\frac{N^{\prime}}{2}\leq\hat{k}\leq N^{\prime}$\end{cases*}$$

$$\begin{array}{rl}{\mathcal{Q}}_1& =\bar{Q}(t_1-t_0)\delta \hat{t}(x_{\hat{k}}-r_{\hat{k}})\\ {\mathcal{Q}}_2& =\bar{Q}(t_f-t_1)\delta \hat{t}(x_{\hat{k}}-r_{\hat{k}})\end{array}$$

$$u^{*}_{\hat{k}}=\begin{cases*}-(\bar{R}(t_{1}-t_{0})\delta\hat{t})^{-1}g_{1}^{T}(x_{\hat{k}})\lambda_{\hat{k}+1}(t_{1},x^{*}_{\hat{k}+1})&if $0\leq\hat{k}<\frac{N^{\prime}}{2}$\\ -(\bar{R}(t_{f}-t_{1})\delta\hat{t})^{-1}g_{2}^{T}(x_{\hat{k}})\lambda_{\hat{k}+1}(t_{1},x^{*}_{\hat{k}+1})&if $\frac{N^{\prime}}{2}\leq\hat{k}\leq N^{\prime}$\end{cases*}$$

$${\hat{\lambda }}_{\hat{k}+1}(t_1,x_{\hat{k}+1})={\hat{W}}_{\hat{k}}^T\varphi (t_1,x_{\hat{k}})$$

$$\widehat{u}_{\hat{k}}(t_{1},x_{\hat{k}})=\begin{cases*}-R_{1}^{-1}g_{1}^{T}\big{(}x_{\hat{k}}\big{)}\widehat{\lambda}_{\hat{k}+1}(t_{1},x_{\hat{k}+1})&if $0\leq\hat{k}<\frac{N^{\prime}}{2}$\\ -R_{2}^{-1}g_{2}^{T}\big{(}x_{\hat{k}}\big{)}\widehat{\lambda}_{\hat{k}+1}(t_{1},x_{\hat{k}+1})&if $\frac{N^{\prime}}{2}\leq\hat{k}\leq N^{\prime}$\end{cases*}$$

Full text

Keywords- Optimal Control, Approximate Dynamic Programming, Autonomous Wheel Loader, Short Loading Cycles.

I Introduction

Wheel loaders (WLs) are essential equipment in off-road constructions. From the control perspective, minimizing fuel consumption and minimizing the operation time of wheel loaders are two crucial and contradictory control goals.

In general, in a WL an engine provides the power. The three actions which consume power in a WL are traction, lifting, and steering. Traction is linked to the driveline for moving the vehicle or stopping it. Lifting is related to the hydraulic system to lift or bring down the boom/bucket, and steering is the one that navigates the vehicle.

In a construction field, wheel loaders mostly go through repetitive cycles to pick up a load from one point and drop it in another point. This cycle includes 4 actions as backwarding, stopping, forwarding and then stopping. After filling the bucket with a load, the vehicle accelerates backward from the resting point. At time $t_{1}$ the operator applies the brakes while putting the gearbox in the neutral mode to bring the vehicle to a full stop at time $t_{2}$. Then the operator puts the gearbox in the forward mode and accelerates toward the unloading point. At time $t_{3}$, the operator applies the brakes again while putting the vehicle in the neutral gearbox to bring the vehicle to full stop at time $t_{4}$. When the vehicle stops, the operator unloads the bucket. The path that a WL follows in such a cycle has a V-shape, and it is called a Short Loading Cycle (SLC).

There are some critical questions in autonomous control of a WL such as

• 1-

What is the best SLC to follow?

• 2-

What are the optimal times to switch from one mode to another mode?

• 3-

What is the optimal time of operation?

• 4-

When is the optimal time to bring the bucket down or lift it?

• 5-

What is the optimal steering angle in an SLC?

• 6-

How can one get the minimum fuel consumption?

In this paper, we tried to address some of these questions in the context of optimal control and left the rest for our future research. The behavior of a WL in an SLC portrays a switched system with controlled subsystems and a fixed mode sequence [1]. Switched systems are systems comprised of several subsystems/modes and a switching rule that assigns the active mode [2, 3, 4]. The sequence of active modes in a switched system is called the mode sequence. For example, consider a system with two modes and only one switching which happens at $t=t_{1}$. One possible mode sequence is $\{mode\hskip 3.0pt1,\hskip 3.0ptmode\hskip 3.0pt2\}$ which means that $\forall t<t_{1}$, mode 1 is active and for $t\geq t_{1}$ mode 2 is active. As seen in this example, when the mode sequence is fixed the controller is only responsible for assigning the switching time instants and not the active mode. Also, when the subsystems include a continuous control, for example steering angle, braking, and hydraulic pressure in a WL, the subsystems are called controlled subsystems. In general optimal control of switched systems is a challenging problem due to discontinuous nature of switching [3, 2, 5, 6].

A complete engine model has many variables/parameters, and it is not suitable for control purpose. In order to model a diesel engine, the data from a real engine [7] is used to tune the parameters of a mean value engine model, introduced in [8], which captures only the essential characteristics of the engine, but maintains a good level of accuracy. After the engine is modeled, the driveline of a WL is modeled as a switching system with forwarding, backwarding, and stopping modes. In the present study, it is assumed that the desired path is known. Hence, an optimal control problem is formulated to track the desired path and to find the optimal switching times.

Compared to the existing literature, the present work neglects the filling and emptying dynamics of the bucket as were studied in [9]. Also, the path is assumed to be known a priori. This is unlike [10, 11, 12, 13, 14] where the path planning is investigated. The present study borrows most of the dynamical modeling from [8]. Also, the current study follows the same idea for modeling the gearbox as a switched system and finding the optimal switching time that was done in [8]. However, the presented optimal control problem formulation and solutions in this paper are entirely different from the ones discussed in [8]. In [8], a mixed integer programming solution provided by a numerical solver is used for a specific initial condition. However, the solution developed in this paper is a closed loop feedback optimal policy which means that the solution is valid for all initial conditions in a compact set selected as a domain of interest. Also, since in this paper the mode sequence is enforced, by a change of mode sequence, one can generate the results reported in [15] in a closed loop feedback policy. Compared to [16], [17], this paper provides the optimal policy in the whole time horizon rather than non-optimal control methods. In [18] optimal control of a wheel loader for gravel application was performed with dynamic programming. In [19], an energy management strategy for hybrid wheel loaders was studied with dynamic programming along with an advisory control policy. In general, dynamic programming is a strong method which provides the closed loop feedback optimal policy [20]. However, as the order of the system increases rapid access to memory becomes prohibitive which is known as the curse of dimensionality [20]. The method provided in this paper is based on Approximate Dynamic Programming (ADP) which finds the near optimal solution instead of the exact optimal solution. This is in fact a remedy for the curse of dimensionality in dynamic programming.

The rest of this paper is organized as follows. In section II, modeling of the diesel engine and the WL are discussed. In section III, optimal control problem formulations are presented first, and then the optimal control solution is discussed. Simulation results are provided in section IV, and section V concludes the paper.

II DYNAMICS OF THE WHEEL LOADER

A model for a wheel loader was introduced in [8] which includes four main parts as powertrain, driveline, steering, and hydraulics. This model can be summarized as follows.

[TABLE]

[TABLE]

In the subsequent subsection, engine model and its components are explained as they were needed. Due to page constraints, the interested readers are referred to [8] for the details of the lifting dynamics and the driveline.

II-A Engine Model

Equations (1) and (2) are the mean value engine model for modeling the engine torque and manifold pressure. In order to find a function for $T_{e}$, the data from a validated diesel engine developed in [7] is used. Hence, random square wave input signals with short ($0.1$ sec) and long ($1.5$ sec) pulse widths were used to excite both transient and steady-state response of the engine. For gathering the training patterns, the engine model was run for $300$ sec subject to short pulse width and $100$ sec with long pulse width. Then the values of $T_{e}$, $\omega_{e}$, $p_{im}$, and $u_{f}$ were first normalized and then saved. For finding a mean value model, $T_{e}:\mathbb{R}^{3}\to\mathbb{R}$ was selected as

[TABLE]

where $W_{T_{e}}\in\mathbb{R}^{4}$ is a tunable weight vector and superscript $T$ denotes the transpose operator. For finding $W_{T_{e}}$ least squares in batch mode was used on the entire data. The Mean Absolute Error (MAE) between the predicted engine torque from (11) and the data measured from the engine for the training was only $7.5\times 10^{-3}$. For validating the model, a new set of random inputs were given to the engine. For gathering the validation data, in the first $50$ sec short pulse width was used and for the second $50$ sec long pulse width was used. The $W_{T_{e}}$ that was found before was used, without retraining, to predict the engine torque. The MAE between the predicted and measured values of engine torque was only $4\times 10^{-2}$ with the validation data which shows the effectiveness of the training with least squares. The performance of the model introduced in (11) and the data from the engine in the validation is shown in Figure 1.

In WLs, the engine provides power for hydraulics, traction, and steering. Hence [8]

[TABLE]

In (12), $P_{lift}$ is the power required for lifting and can be defined as

[TABLE]

In (13), $Q(.)$ is the mass flow rate of the hydraulic fluid into the lifting cylinder and $\eta_{lift}$ is the lifting efficiency. Also, steering power can be shown as

[TABLE]

where $C_{st}$ is a constant. For modeling the traction power, [8] proposed a hybrid model based on the gear ratio as

[TABLE]

Assuming three values as $\gamma\in\{-60,0,+60\}$, one can consider the switching dynamics as backwarding, stopping, and forwarding.

For identifying the manifold pressure, the same approach as using the data from the validated engine model was used to tune the parameters of the model introduced in [21], [22] as

[TABLE]

[TABLE]

In (16) and (17), $\tau_{1}$, $\tau_{2}$, $a_{1}$, $a_{2}$, and $a_{3}$ are constants. For this tuning, gradient descent algorithm was used with the sequential offline data. Once the training concluded, the model was validated with a new set of data. The MAE between the trained model and the validation data was only $5.3\times 10^{-2}$.

II-B State space model

It is desired to show the dynamics of the WL as

[TABLE]

where $x\in\mathbb{R}^{n}$ is the state vector, and Lipschitz functions $\bar{f}_{v}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ and $\bar{g}_{v}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ denote the dynamics of the subsystems. The active mode in time instant $t$ is shown by sub-index $v$, and the set of all subsystems is shown as $\{1,2,3\}$ which corresponds to forwarding, backwarding, and stopping mode in the WL. Equation (18) shows a switched system dynamics where each subsystem includes a continuous control. In the literature of the switched systems, this type of subsystems is called controlled subsystems. Also, equation (18) shows that each mode is control affine. To show the dynamics of the WL as (18), let the state vector as

[TABLE]

where time dependency of $x_{1\leq i\leq 9}$ is dropped for notational simplicity, i.e., $x_{i}=x_{i}(t)$. Also, $x_{1}=\omega_{e}$, $x_{2}=p_{im}$, $x_{3}=\theta$, $x_{4}=\omega$, $x_{5}=X$, $x_{6}=Y$, $x_{7}=V$, $x_{8}=\beta$, $x_{9}=\delta$. In [8], the control vector was selected as $u(t)=[u_{f},u_{p},u_{b},u_{s}]^{T}$. In this research we are interested to find the control affine dynamics and existence of nonlinear terms with respect to $u_{s}$ and $u_{p}$ will not permit such selection for $u(t)$. To solve this problem, similar to the method used in [23], two new states can be defined as $x_{10}=u_{p}$ and $x_{11}=u_{s}$. Taking time derivative from $x_{10}$, and $x_{11}$, one has

[TABLE]

Also, one can select $u_{3}=u_{f}$ and $u_{4}=u_{b}$. Therefore, $u(t)=[u_{1},u_{2},u_{3},u_{4}]^{T}$.

Remark 1

In order to normalize the range of variations of the variables, one can nondimensionalize the dynamics. For this purpose, let $X_{1}$ as the maximum value for variable $|x_{1}|$. Hence, one can define $\bar{x}_{1}=\frac{x_{1}}{X_{1}}$. With this transformation, one can see that $|x_{1}|\in[0,X_{1}]$ and $|\bar{x}_{1}|\in[0,1]$. Considering a scalar system as $\dot{x_{1}}=f(x_{1})+g(x_{1})u(t)$, it is straight forward to see \dot{\bar{x}}_{1}=\frac{1}{X_{1}}\big{(}f(X_{1}\bar{x}_{1})+g(X_{1}\bar{x}_{1})u(t)\big{)}. The discussed nondimensionalization is used in this paper for optimal controller design.

III OPTIMAL CONTROL FORMUATION

Considering the dynamics as (18) and assuming the mode sequence is known, it is desired to find optimal switching times and a feedback control policy that minimizes the cost function

[TABLE]

In (22), $t_{0}$ is the initial time, $t_{f}$ is the final time, and $r\in\mathbb{R}^{n}$ is the reference signal which is a known function of time. $S\in\mathbb{R}^{n\times n}$ is a positive semi-definite matrix for penalizing the terminal cost, $\bar{Q}\in\mathbb{R}^{n\times n}$ is the state penalizing matrix which is assumed to be positive semi-definite, and $\bar{R}\in\mathbb{R}^{m\times m}$ is a positive definite control penalizing matrix.

III-A Including Mode Sequence

To solve the optimal switching problem, one needs to solve a two-level optimization [24]. In the upper level, switching times are assigned and in the lower level, the control policies to ensure the tracking are sought. To introduce the idea, consider a switched system with two subsystems and only one switching which happens at $t=t_{1}$. Also, let the mode sequence be $\{mode$ $1,$ $mode$ $2\}$. To make the switching time instant an independent parameter, let [24]

[TABLE]

using chain rule one can find $x^{\prime}(\hat{t})=\frac{dx}{d\hat{t}}=\frac{dx}{dt}\frac{dt}{d\hat{t}}$ as

[TABLE]

Also, changing the independent variable in the cost function (22) from $t$ to $\hat{t}$ leads

[TABLE]

Letting $\delta\hat{t}$ as a small sampling time, by using Euler method one can discretize (24) and (25) as

[TABLE]

[TABLE]

In (26), f_{1}\big{(}x_{\hat{k}}\big{)}=x_{\hat{k}}+\bar{f}_{1}\big{(}x_{\hat{k}}\big{)}(t_{1}-t_{0})\delta\hat{t}, g_{1}\big{(}x_{\hat{k}}\big{)}=\bar{g}_{1}\big{(}x_{\hat{k}}\big{)}(t_{1}-t_{0})\delta\hat{t}, f_{2}\big{(}x_{\hat{k}}\big{)}=x_{\hat{k}}+\bar{f}_{2}\big{(}x_{\hat{k}}\big{)}(t_{f}-t_{1})\delta\hat{t}, and g_{2}\big{(}x_{\hat{k}}\big{)}=\bar{g}_{2}\big{(}x_{\hat{k}}\big{)}(t_{f}-t_{1})\delta\hat{t}. Also, in (26), $\hat{k}\in[1,N^{\prime}]$ is the discrete time index where $N^{\prime}=\frac{number\hskip 3.0ptof\hskip 3.0ptswitching+1}{\delta\hat{t}}$ [1]. For finding the minimum cost-to-go from discrete time index $\hat{k}$ to $N^{\prime}$(value function) one has111For $\hat{k}<N^{\prime}/2$. Otherwise, $(t_{1}-t_{0})$ in (28) should be replaced by $(t_{f}-t_{1})$.

[TABLE]

The optimal policy can be defined as

[TABLE]

Taking the gradient of the optimal value function, one can define the optimal costate as

[TABLE]

where

[TABLE]

It is straight forward to see if optimal costates are known, one can directly calculate the optimal policy as

[TABLE]

In (30) and (32), the super-script $*$ denotes the optimality. In the next subsection, a solution is presented to find the optimal costates.

III-B ADP Solution

The backbone of the solution is training neural networks to approximate $\lambda_{\hat{k}+1}(t_{1},x_{\hat{k}+1})$ from $(t_{1},x_{\hat{k}})$. Based on Weierstrass Approximation Theorem [25], linear-in-parameter neural networks with polynomial basis functions can uniformly [26] approximate continuous functions to a desired degree of precision in a compact set. Assuming that the value functions are continuously differentiable, one can use linear-in-parameter neural networks to approximate the costates. Hence, the approximate costates can be calculated as

[TABLE]

where $\widehat{W}_{\hat{k}}\in\mathbb{R}^{m_{\lambda}}\times\mathbb{R}^{n}$ is a tunable weight vector and $m_{\lambda}$ is the number of polynomial basis functions (neurons). The weight vectors are tuned through the training process backward in time. Once the costates are known, one finds the approximate optimal policy as

[TABLE]

In (34), $R_{1}=\bar{R}\delta\hat{t}(t_{1}-t_{0})$ and $R_{2}=\bar{R}\delta\hat{t}(t_{f}-t_{1})$. As mentioned before, for training, one can go backward in time and find the costates and save them for online control.

Once the training concluded, one needs to find the optimal switching times from the costates for a selected initial condition $x_{0}\in\Omega$. For this purpose, one can propagate the states along all possible switching times by using the costates and find the optimal cost to go for all possible switching time. Once done, one can choose the switching times which lead to the minimum cost.

IV SIMULATION RESULTS

To start the simulations, the values selected for the parameters of the WL model are given in Table II. Some of these parameters are not introduced in this paper due to page constraints. Interested readers are directed to [8] for more details.

Also, $W_{T_{e}}=[-0.154712845456646$, $0.555616949283938,$ $4.84689263976440,$ $0.450411638112411]^{T}$ in (11) and the parameters in (16) and (17) are given as $a_{1}=3.07399910699804$, $a_{2}=0.467292644030092$, $a_{3}=0.467292644032403$, $\tau_{1}=5.94705195850280$, $\tau_{2}=-0.703974251212264$.

In order to see the performance of the model, some selected control inputs are given to the WL, and the performance of the WL is evaluated. The mode sequence in this example is $\{$backwarding, stopping, forwarding, stopping$\}$. All the switching time instants and the continuous control inputs to the wheel loader are selected manually to simulate a short loading cycle. The braking input, vehicle velocity and the path are shown in Figures 2, 3, and 4, respectively.

In order to apply the ADP method to find the optimal switching time and control, one assumes that the reference signal is a known function of time. For this example, the reference signal is selected as $\dot{r}(t)=[0,0,0,0,\sin(\pi t),\sin(\pi t),0,0,0,0,0]$. The start and end times are selected as $t_{0}=0$, and $t_{f}=3$. In this example, the mode sequence is selected as $\{backwarding,\hskip 3.0ptstopping\}$. In other words, half of an SLC is considered. Since there is only one switching in the system, by choosing $\delta\hat{t}=0.001$, one can calculate $N^{\prime}=2000$. The state and control penalizing matrices are selected as $S=\operatorname{diag}([0,0,0,0,10^{4},10^{4},10^{4},0,0,0,0])$, $\bar{Q}=\operatorname{diag}([0,0,0,0,10^{4},10^{4},0,0,0,0,0])$, and $\bar{R}=\operatorname{diag}([1,1,1,1])/\delta\hat{t}$ where $\operatorname{diag}([a,b])$ is a diagonal matrix with values of $a$ and $b$ on the main diagonal and zero elsewhere. In choosing $S$ matrix for terminal state penalizing, one notes that penalization of the velocity is enforced as it is important to bring the vehicle to rest at the final time.

In order to start the training, one needs to select a good neural network to capture the dynamics of the costates. Mostly, such selection is performed by trial and error. After selecting a relatively rich set of basis functions, the training process was performed. The history of the neural network weights is shown in Figure 5. As one can see from Figure 5, the history of the weights shows a jump at $\hat{t}=1$ which is the switching time.

When the training process concluded, a random initial condition was selected in the domain of training. The initial condition for the reference signal was selected as $r_{0}=[0,0,0,0,0.3,-0.1,0,0,0,0,0]^{T}$ and the optimal switching time was found as $t_{1}=2.86$. Afterward, the optimal costates were used to propagate the states from the previously selected random initial condition. The history of the states is shown in Figure 6. The poor tracking is the result of approximation errors in capturing the dynamics of the costates and it is ongoing research with our research group to solve this problem.

V CONCLUSION

An approximate dynamic programming solution is used for optimal control of a wheel loader. For this purpose, the dynamics of the wheel loader are modeled as a switched system with controlled subsystems and a fixed mode sequence. Some simulation results are provided to show the effectiveness of the solution.