Published: 01 August 2023

# Minimum-time lane changing problem of vehicle handling inverse dynamics based on adaptive mesh refinement and collocation optimization method

Yingjie Liu1
Dawei Cui2
Wen Peng3
1, 2School of Machinery and Automation, Weifang University, Weifang, 261061, Shandong, China
3State Key Laboratory of Rolling and Automation, Northeastern University, Shenyang, 110819, China
Corresponding Authors:
Yingjie Liu, Dawei Cui
Views 26

#### Abstract

In order to solve the problems of multiple constraints, many different calculations of nonlinear equations which lead to major errors in the process of vehicle lane changing with minimum time, an adaptive mesh refinement and collocation optimization method is proposed. Firstly, the problem of vehicle lane changing with minimum time has been divided into nonlinear programming problems in different grids. The Lagrange interpolation polynomial was used to approximate the solution of the optimization problem in the grid, and the absolute and relative errors were resolved. Then, a mesh was determined for the unsmooth part according to the curvature of the trajectory, and the location and number of meshes were further determined according to the relationship between the maximum relative error and the allowable error. At the same time, the solution accuracy was improved by adding an adaptive calculation to the smooth interval which does not meet the tolerance error. Finally, the simulation example of comparison with the traditional optimization method was proposed. The results showed that the algorithm presented in the paper had a higher solution efficiency under the same calculation accuracy. #### Highlights

• The method can adaptively increase the number of collocation and reduce the number of mesh.
• The method can generate fewer total collocation in the solution area reducing the number of nonlinear programming equations after discretization.
• The algorithm presented in the paper has a higher solution efficiency under the same calculation accuracy.

## 1. Introduction

The rapid development of automobile technology has brought great convenience to human travel and life, and has also driven the rapid growth of automobile ownership in recent years. At present, automobile has become one of the pillar industries of national economy in various countries. However, due to some improper operations of drivers in complex traffic scenes, more and more traffic accidents seriously threaten the safety of human life. The longitudinal and lateral movements of vehicles are not completely independent, but there is a nonlinear coupling relationship. In some cases, the longitudinal and lateral coupling nonlinear characteristics will be very significant, and will affect the accuracy of vehicle motion control and even tracking, and even cannot ensure the vehicle stability in serious cases. For example, when an intelligent vehicle encounters sudden dangerous situations such as sudden braking in front of a static obstacle in the same lane or sudden emergency braking in front of the preceding vehicle, it can learn from an experienced driver to make a lateral lane change with a large range of steering operations. And at the same time, a certain braking operation is introduced to better increase the distance from the dangerous obstacles. At this time, the motion state of the intelligent vehicle traveling with large longitudinal and lateral accelerations may exceed the boundary of the power domain corresponding to the steady state condition, leading the vehicle to enter a strong transient condition with significant nonlinear dynamic characteristics of longitudinal and lateral coupling. Under such a strong transient condition, the longitudinal component of the lateral force of the vehicle will aggravate the impact on the longitudinal speed of the vehicle. And the change of the longitudinal speed will also cause the change of the centrifugal force and the vertical load, causing interference to the lateral motion state of the vehicle. With the action range of longitudinal and lateral acceleration further approaching to the limit boundary value of the dynamic domain provided by road adhesion, the vehicle will obtain a critical instability condition. In order to ensure the control margin of the vehicle, the intelligent driving system usually limits the lateral acceleration in its trajectory plan. Therefore, under the critical instability condition, the vehicle is often accompanied by a large longitudinal acceleration, which causes a large change in the longitudinal speed, leading to a more significant nonlinear dynamics of the vehicle longitudinal and lateral coupling. As the current control strategy considers the longitudinal and lateral designs of intelligent vehicle motion separately and lacks comprehensive consideration of the above longitudinal and lateral coupling nonlinear dynamic characteristics, the vehicle motion control performance may decline significantly under strong transient, critical instability and other working conditions, which lead to vehicle instability too .

In the field of intelligent vehicle research, the content related to autonomous lane changing has always been one of the research hotspots in the industry. According to incomplete statistics, vehicle traffic accidents when vehicles change lanes occur every year in the quantity of about 15 %-20 %. In addition, traffic jams caused by changing lanes account for about 25 % of the total. Lane changing is a complex vehicle behavior, which not only affects the safety, efficiency and comfort of the vehicle itself, but also greatly affects the performance of the entire traffic system as a part of the traffic flow. Therefore, it is of great significance to develop a set of intelligent vehicle automatic lane changing system to meet multi-objective optimization. From the current market situation, a very mature vehicle automatic lane change system that has reached the mass production level has not yet appeared .

The problem of vehicle lane change motion planning has been widely studied. A brief review is presented in the following.

Naskath et al.  analyzed the connectivity of high-speed mobility and lane changing based on discretionary lane changing approach for V2V environment. Muhamad et al.  proposed a yaw rejection control for a single-trailer truck using a steerable wheel located at the middle axle of the truck. Tolga et al.  analyzed autonomous electric vehicles and compared their potential environmental impacts with other public transport types, carpooling, walking, cycling, and various transportation policy applications. Cao et al.  designed a lane change algorithm based on a look-ahead concept for the vehicle driving in front of the emergency vehicle to avoid emergency situations. Jiang et al.  presented a lane-level vehicle counting system which was based on V2X communications and centimeter-level positioning technologies. Yang et al.  presented a road environment representation method based on a dynamic occupancy grid for the lane changing assistance strategy. Vaishnavi et al.  developed a lane optimization model for the CAV system for solving the complexities of multilane merging areas. Amit et al.  proposed a microscopic traffic flow model which was intended to explain accurately the lane changing activity using Cellular Automata. Xiao et al.  proposed a non-lane-discipline-based car-following model for solving the problem of vehicle lane changing. Huang et al.  proposed a neural network-based method to predict the future status of the local vehicle using the information from V2X, and estimate the future energy consumption of each lane. Panagiotis et al.  proposed a methodology for modeling the likelihood of lane changing at intersections to quantify the favorability of the surrounding environment towards lane changing and simple LSTM modeling structures. Wang et al.  adopted a generalized dynamic model for automatic lane-changing to solve the limitations of macroscopic evaluation. H. Hamedi et al.  studied the vehicles’ lane change trajectories to present three novel prediction models for trajectory-grounded position. In order to investigate the characteristics of double fire sources, Guo et al.  performed four full-scale fire experiments in the Jiangpu road tunnel. Sharma et al.  proposed a novel hierarchical software architecture for the prediction of lane changing behavior on highways. Zhou et al.  proposed a hierarchical control strategy for the vehicle stability control. Wang et al.  proposed a two-stage crash causation analysis method based on the pre-crash scenarios and a crash causation derivation framework that systematically categorized and analyzed contributing factors. Wei et al.  proposed a prediction model based on an attention-aided encoder-decoder structure and deep neural network (DNN) to solve the problems of low prediction accuracy, difficulty in long-term prediction. Zhang et al.  carried out a naturalistic driving study on a highway, from which a large amount of on-road data consisting of lane-keeping and lane-changing left and right maneuvers were collected. Tajalli et al.  presented a methodology for optimal control of connected automated vehicles in freeway segments with a lane drop. He et al.  proposed a novel observation adversarial reinforcement learning approach for robust lane change decision making of autonomous vehicles. Shi, et al.  proposed an integrated deep learning-based two-dimension trajectory prediction model which could predict combined behaviors. Li et al.  proposed a lane change decision-making framework based on deep reinforcement learning to find a risk-aware driving decision strategy with the minimum expected risk for autonomous driving to ensure the driving safety.

Some of the methods provided in the above documents have large solution errors, some cannot converge to the optimal solution, and some have low accuracy. In order to solve the above problems, the relative demerits were first eliminated in discrete grids based on the hp-adaptive pseudospectral method as a direct method. Then, the mesh to be refined is determined according to the curvature of state and control variables. Finally, according to the relationship between relative and tolerance error, the number of sub grids and the number of collocations in the grid are determined, and then a solution method for vehicle lane changing with minimum time is constructed. A typical example is used to verify the effectiveness of the method. In the process of solving minimum time lane changing problem of vehicle handling inverse dynamics by a direct method, a large number of nodes need to be added to improve the accuracy of solution. So, it will lead to a dense Jacobian matrix for solving nonlinear programming problems, high order of interpolation polynomial, large amount of required calculations. Moreover, the increase of nodes will also lead to the poor convergence of the algorithm. Therefore, it is necessary to optimize the mesh and collocation. For the problem of minimum time lane changing of vehicle handling inverse dynamics, the number of sub-grid divisions is provided by using the relative error, and the location of sub-grid divisions is determined by using the density function. The adaptive mesh and refinement and collocation optimization method can approach the optimal solution by densifying the mesh and adding collocation points nearby. In areas with high smoothness, fewer collocations can be used with the strong adaptability. In the process of optimization, the mesh should not be refined for the smooth parts that do not meet the requirements of allowable error, while the collocation should be added. The method can adaptively increase the number of collocation and reduce the number of meshes, so that the total number of collocations generated in the solution area is less, reducing the number of discrete nonlinear programming equations, improving the speed of solving the optimal trajectory by the SQP method.

## 2. Mathematical model of vehicle lane change problem

### 2.1. Mathematical model of vehicle lane change problem

A nonlinear 7-DOF vehicle model shown in Fig. 1 is used to describe the vehicle lane change problem.

In the state space form, it is :

1
$\left\{\begin{array}{l}m\left({\stackrel{˙}{v}}_{x}-r{v}_{y}\right)=\left({F}_{xfl}+{F}_{xfr}\right)\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_{f}+{F}_{xrl}+{F}_{xrr}-\left({F}_{yfl}+{F}_{yfr}\right)\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_{f}\right),\\ m\left({\stackrel{˙}{v}}_{y}+r{v}_{x}\right)=\left({F}_{xfl}+{F}_{xfr}\right)\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_{f}+{F}_{yrl}+{F}_{yrr}+\left({F}_{yfl}+{F}_{yfr}\right)\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_{f}\right),\\ {I}_{z}\stackrel{˙}{r}=\left({F}_{xfl}+{F}_{xfr}\right){l}_{f}\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_{f}+\left[\left({F}_{xfr}-{F}_{xfl}\right)\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_{f}+\left({F}_{xrr}-{F}_{xrl}\right)\right]\frac{{B}_{w}}{2}\\ +\left({F}_{yfl}+{F}_{yfr}\right){l}_{f}\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_{f}+\left({F}_{yfl}-{F}_{yfr}\right)\frac{{B}_{w}}{2}\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_{f}-\left({F}_{yrl}+{F}_{yrr}\right){l}_{f}.\end{array}\right\$

The state equation can be described as:

2
$\stackrel{˙}{x}=f\left[x\left(t\right),z\left(t\right)\right],$

where $x\left(t\right)$ and $z\left(t\right)$ are the state and the input which are denoted respectively as $x\left(t\right)={\left[{v}_{x}\left(t\right),{v}_{y}\left(t\right),r\left(t\right),x\left(t\right),y\left(t\right)\right]}^{T}$, $z\left(t\right)=\left[\delta \left(t\right){\right]}^{T}$.

The parameters and the corresponding definitions can be found in Table 1.

Fig. 1Nonlinear 7-DOF vehicle model Table 1Parameter and definition

 Parameter Definition ${v}_{x}$ Longitudinal speed ${v}_{y}$ Lateral speed $r$ Yaw rate of vehicle $m$ Vehicle mass ${I}_{z}$ Moment of inertia around z axis $\beta$ Sideslip angle ${B}_{w}$ Wheel-base of vehicle ${l}_{f}$ Distances of front axle from center of gravity ${l}_{r}$ Distances of rear axle from tcenter of gravity ${F}_{xij}$ Longitudinal tire force ${F}_{yij}$ Lateral tire force ${\delta }_{f}$ Front steering angle where $i=f$ or $r$, $j=l$ or $r$, and $fl$ means the front left, $fr$ means the front right, $rl$ means the rear left, $rr$ means the rear right

### 2.2. Tire model

3
$y=D\mathrm{s}\mathrm{i}\mathrm{n}\left\{C\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left[Bx-E\left(Bx-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(Bx\right)\right)\right]\right\}+S,$
4
$Y\left(X\right)=y\left(x\right)+{S}_{v},$
5
$x=X+{S}_{h}.$

The parameters and the corresponding definitions in Eqs. (3)-(5) can be found in Table 2.

Table 2Parameter and definition

 Parameter Definition $B$ Stiffness factor $C$ Shape factor $D$ Amplitude factor $E$ Curvature factor ${S}_{h}$ Curve offset on the $x$-axis ${S}_{v}$ Curve offset on the $y$-axis $X$ Input variable $Y\left(X\right)$ Output variable

### 2.3. Optimal control object of lane change motion planning problem

The minimum time problem of lane change motion planning can be regarded as the optimal control problem in the control theory. Therefore, the minimum time performance indicator is:

6
$J={\int }_{{t}_{0}}^{{t}_{f}}dt,$

where:

${D}_{ki}\left({\tau }_{k}\right)=\stackrel{˙}{L}\left({\tau }_{k}\right)=\left\{\begin{array}{l}\frac{\left(1+{\tau }_{k}\right){\stackrel{˙}{P}}_{N}\left({\tau }_{k}\right)+{P}_{N}\left({\tau }_{k}\right)}{\left({\tau }_{k}-{\tau }_{i}\right)\left[\left(1+{\tau }_{i}\right){\stackrel{˙}{P}}_{N}\left({\tau }_{i}\right)+{P}_{N}\left({\tau }_{i}\right)\right]},i\ne k,\\ \frac{\left(1+{\tau }_{i}\right){\stackrel{¨}{P}}_{N}\left({\tau }_{i}\right)+2{\stackrel{˙}{P}}_{N}\left({\tau }_{i}\right)}{2\left[\left(1+{\tau }_{i}\right){\stackrel{˙}{P}}_{N}\left({\tau }_{i}\right)+{P}_{N}\left({\tau }_{i}\right)\right]},i=k,\end{array}\right\$

is the initial time, ${D}_{ki}\in {R}^{N×\left(N+1\right)}$ is the final time.

### 2.4. Constrains

The initial and terminal states are described as:

7
$x\left({t}_{0}\right)=\left[{v}_{x}\left({t}_{0}\right),0,0,0,0{\right]}^{T},$
8
$x\left({t}_{f}\right)=\left[{v}_{x}\left({t}_{f}\right),0,0,x\left({t}_{f}\right),y\left({t}_{f}\right){\right]}^{T}.$

In order to accomplish the lane change motion planning maneuver successfully, the constraints set on the longitudinal and lateral lines as well as the lateral acceleration are $0\le x\le {x}_{b}$, ${B}_{0}\le y\le {B}_{1}$, $\left|{a}_{y}\right|\le$ 3 m/s2. Where ${x}_{b}$ is the longitudinal distance between vehicle and obstacle; ${B}_{0}$and${B}_{1}$ are the minimum lateral displacements required for completing lane change motion planning process and the left boundary of adjacent lane respectively .

When the braking maneuver is applied to decelerate the vehicle, the constraints on ${F}_{xf}$, ${F}_{xr}$ can be rewritten in the following manner:

9
$\left\{\begin{array}{l}{F}_{xf}\ge -\frac{\mu mg\left(b+\mu {h}_{g}\right)}{a+b},\\ {F}_{xr}=\frac{a-\mu {h}_{g}}{b+\mu {h}_{g}}{F}_{xf}.\end{array}\right\$

## 3. Lane changing with minimum time

### 3.1. Problem description

For the sake of convenience, the problem of Lane Changing with Minimum Time is transformed into a Mayer problem.

The Bolza cost function is:

10
$\begin{array}{l}J=\psi \left(x\left(-1\right),{t}_{0},x\left(1\right),{t}_{f}\right)+\frac{{t}_{f}-{t}_{0}}{2}{\int }_{-1}^{1}g\left(x\left(\tau \right),z\left(\tau \right),\tau ;{t}_{0},{t}_{f}\right)d\tau ,\\ s.t.\left\{\begin{array}{l}\frac{dx}{d\tau }=\frac{{t}_{f}-{t}_{0}}{2}f\left(x\left(\tau \right),z\left(\tau \right),\tau ;{t}_{0},{t}_{f}\right),\tau \in \left[-1,1\right],\\ E\left(x\left(-1\right),{t}_{0},x\left(1\right),{t}_{f}\right)=0,\\ C\left(x\left(\tau \right),z\left(\tau \right),\tau ;{t}_{0},{t}_{f}\right)\le 0,\tau \in \left[-1,1\right],\end{array}\right\\end{array}$

where $J$ is the cost function; $\psi \left(\cdot \right)$ is the Mayer type cost function; $g\left(\cdot \right)$ is the Lagrange type cost function; $f\left(\cdot \right)$ indicates the dynamic equation; $E\left(\cdot \right)$ is the boundary constraint; $C\left(\cdot \right)$ is the path constraint; $x\in {R}^{n}$ and $z\in {R}^{m}$ represent the state and control variables respectively; ${t}_{0}$ is the initial time; ${t}_{f}$ is the final time; $x\left(-1\right)$ and $x\left(1\right)$ are the initial and final values of the state variable. Because the solution interval is [–1, 1], $t$ should be transformed as follows:

11
$t=\tau \left({t}_{f}-{t}_{0}\right)/2+\left({t}_{f}+{t}_{0}\right)/2.$

In the discrete process using the HP method, the solution interval $\tau \in \left[-1,1\right]$ will be divided into $N$ grids, i.e. ${S}_{k}=\left[{T}_{k-1},{T}_{k}\right]$, $k=1,2,\cdots ,N$, meeting $-1={T}_{0}<{T}_{1}<\cdots <{T}_{k}=1\text{,}$${\bigcup }_{k=1}^{N}{S}_{k}=\left[-1,1\right]$.

It is set that ${x}^{\left(k\right)}\left(\tau \right)$ and ${z}^{\left(k\right)}\left(\tau \right)$ represent the state and control variables of ${S}_{k}$ respectively. Then the above optimal control Bolza problem can be transformed as:

12
$\begin{array}{l}J=\psi \left({x}^{\left(1\right)}\left(-1\right),{t}_{0},{x}^{N}\left(1\right),{t}_{f}\right)+\frac{{t}_{f}-{t}_{0}}{2}\sum _{k=1}^{N}{\int }_{{T}_{k-1}}^{{T}_{k}}g\left({x}^{\left(k\right)}\left(\tau \right),{z}^{\left(k\right)}\left(\tau \right),\tau ;{t}_{0},{t}_{f}\right)d\tau ,\\ s.t.\left\{\begin{array}{l}\frac{d{x}^{\left(k\right)}\left(\tau \right)}{d\tau }=\frac{{t}_{f}-{t}_{0}}{2}f\left({x}^{\left(k\right)}\left(\tau \right),{z}^{\left(k\right)}\left(\tau \right),\tau ;{t}_{0},{t}_{f}\right),k=1,2,\cdots ,N,\\ E\left({x}^{\left(1\right)}\left(-1\right),{t}_{0},{x}^{\left(N\right)}\left(1\right),{t}_{f}\right)=0,\\ C\left({x}^{\left(k\right)}\left(\tau \right),{z}^{\left(k\right)}\left(\tau \right),\tau ;{t}_{0},{t}_{f}\right)\le 0,k=1,2,\cdots ,N.\end{array}\right\\end{array}$

Because the state is continuous at the end of the grid, $x\left({T}_{k-1}^{+}\right)=x\left({T}_{k}^{-}\right)$ should be met at the connection point.

The pseudo-spectral method can be divided into four categories according to different orthogonal polynomials. But the discrete principle is basically the same among all polynominals. The paper takes the Radau pseudo-spectral method as an example to illustrate discrete principle of the pseudo-spectral method.

The Bolza optimal control problem is discretized on ${S}_{k}$ grid using LGR (Legendre-Gauss-Radau) collocation points :

13
${x}^{\left(k\right)}\left(\tau \right)\approx {X}^{\left(k\right)}\left(\tau \right)=\sum _{j=1}^{{M}_{k}+1}{X}_{j}^{\left(k\right)}\left({\tau }_{i}\right){L}_{j}^{\left(k\right)}\left(\tau \right),$

where:

14
${L}_{j}^{\left(k\right)}\left(\tau \right)=\prod _{\begin{array}{l}l=1\\ l\ne j\end{array}}^{{M}_{k}+1}\frac{\tau -{\tau }_{l}^{\left(k\right)}}{{\tau }_{j}^{\left(k\right)}-{\tau }_{l}^{\left(k\right)}},\tau \in \left[-1,+1\right],$

where ${L}_{j}^{\left(k\right)}\left(\tau \right)$$\left(j=1,2,\cdot \cdot \cdot ,{M}_{k}+1\right)$ is the Lagrange interpolating polynomials; $\left({\tau }_{1}^{\left(k\right)},{\tau }_{2}^{\left(k\right)},\cdots ,{\tau }_{{M}_{k}}^{\left(k\right)}\right)$ are the collocation points of LGR in interval ${S}_{k}=\left[{T}_{k-1},{T}_{k}\right]$. But ${\tau }_{{M}_{k}}^{\left(k\right)}+1={T}_{k}$ is not the collocation point. ${X}^{\left(k\right)}$ is the approximation of ${x}^{\left(k\right)}$.

Eq. (15) can be obtained by operating difference to $\tau$:

15
$\frac{d{X}^{k}\left(\tau \right)}{d\tau }=\sum _{j=1}^{{M}_{k}+1}{X}_{j}^{\left(k\right)}\frac{d{L}_{j}^{\left(k\right)}\left(\tau \right)}{d\tau }.$

Approximation of the dynamic equation on the grid ${S}_{k}$ is:

16
$\sum _{j=1}^{{M}_{k}+1}{D}_{ij}^{\left(k\right)}{X}_{j}^{\left(k\right)}-\frac{{t}_{0}-{t}_{f}}{2}f\left({X}_{i}^{\left(k\right)}\left(\tau \right),{Z}_{i}^{\left(k\right)}\left(\tau \right),t\left({\tau }_{i}^{\left(k\right)},{t}_{0},{t}_{f}\right)\right)=0,$

where:

17
${D}_{ij}^{\left(k\right)}=\frac{d{l}_{j}^{\left(k\right)}\left({\tau }_{i}\right)}{d\tau },i=1,2,\cdots ,{M}_{k},j=1,2,\cdots ,{M}_{k}+1,$

where ${D}_{ij}^{\left(k\right)}$ is the ${M}_{k}×\left({M}_{k}+1\right)$-dimensional LGR differential matrix on grid ${S}_{k}$. Then Eq. (12) can be discretized as:

18
$\begin{array}{l}J=\psi \left({X}_{\left(1\right)}^{\left(1\right)},{t}_{0},{X}_{{}^{\left({M}_{k}+1\right)}}^{\left(1\right)},{t}_{f}\right)+\sum _{k=1}^{N}\sum _{j=1}^{{M}_{k}}\frac{{t}_{k}-{t}_{k-1}}{2}{\omega }_{j}^{\left(k\right)}g\left({X}_{j}^{\left(k\right)}\left(\tau \right),{Z}_{j}^{\left(k\right)}\left(\tau \right),t\left({\tau }_{j}^{\left(k\right)},{t}_{0},{t}_{f}\right)\right),\\ s.t.\left\{\begin{array}{l}\sum _{j=1}^{{M}_{k}+1}{D}_{ij}^{\left(k\right)}{X}_{j}^{\left(k\right)}-\frac{{t}_{k}-{t}_{k-1}}{2}f\left({X}_{i}^{\left(k\right)}\left(\tau \right),{Z}_{i}^{\left(k\right)}\left(\tau \right),t\left({\tau }_{i}^{\left(k\right)},{t}_{0},{t}_{f}\right)\right)=0,\\ E\left({X}^{\left(1\right)}\left(-1\right),{t}_{0},{X}^{\left(N\right)}\left(+1\right),{t}_{N}\right)=0,\\ \frac{{t}_{k}-{t}_{k-1}}{2}C\left({X}_{i}^{\left(k\right)}\left(\tau \right),{Z}_{i}^{\left(k\right)}\left(\tau \right),t\left({\tau }_{i}^{\left(k\right)},{t}_{k},{t}_{k-1}\right)\right)\le 0,\\ {X}_{{M}_{k}+1}^{\left(k-1\right)}={X}_{1}^{\left(k\right)},\end{array}\right\\end{array}$

where $i=1,2,\cdots ,{M}_{k}$; $k=1,2,\cdots ,N$.

According to the above transformation method, the optimal control problem is transformed into a constrained nonlinear programming problem (NLP), which can be solved by the Sequential Quadratic Programming (SQP) method.

## 4. Adaptive optimization of mesh refinement and collocation

When the smoothness of the state and control variables is low, the solution using the global interpolation polynomial will affect its accuracy greatly. In order to improve the efficiency and accuracy of the solution, it is necessary to solve the discrete error to determine whether it is necessary or not to mesh the interpolated nodes.

### 4.1. Discrete error solution

It is assumed that the nonlinear programming problem shown in Eq. (10) has been already solved in interval ${S}_{k}=\left[{T}_{k-1},{T}_{k}\right]$, $k\in \left[1,2,\cdots ,N\right]$. There are ${M}_{k}+1$ LGR collocation points i.e. $\left({\stackrel{^}{\tau }}_{1}^{\left(k\right)},\cdots ,{\stackrel{^}{\tau }}_{{M}_{k}}^{\left(k\right)},{\stackrel{^}{\tau }}_{{M}_{k}+1}^{\left(k\right)}\right)$ in ${S}_{k}=\left[{T}_{k-1},{T}_{k}\right]$, where ${\stackrel{^}{\tau }}_{1}^{\left(k\right)}={\tau }_{1}^{\left(k\right)}={T}_{k-1}$, ${\stackrel{^}{\tau }}_{{M}_{k}+1}^{\left(k\right)}={T}_{k}$. It is set that the approximation of state variables in $\left({\stackrel{^}{\tau }}_{1}^{\left(k\right)},\cdots ,{\stackrel{^}{\tau }}_{{M}_{k}}^{\left(k\right)},{\stackrel{^}{\tau }}_{{M}_{k}+1}^{\left(k\right)}\right)$ is $\left(X\left({\stackrel{^}{\tau }}_{1}^{\left(k\right)}\right),\cdots ,X\left({\stackrel{^}{\tau }}_{{M}_{k}}^{\left(k\right)}\right),X\left({\stackrel{^}{\tau }}_{{M}_{k}+1}^{\left(k\right)}\right)\right)$. $Z\left({\stackrel{^}{\tau }}_{1}^{\left(k\right)}\right)$ is the approximation of control variables in ${\tau }_{i}^{\left(k\right)}$.

Additionally, the control in ${S}_{k}$ is approximated using Lagrange interpolating polynomials as:

19
${Z}^{\left(k\right)}\left(\tau \right)=\sum _{j=1}^{N}{Z}_{j}^{\left(k\right)}\left(\tau \right){\stackrel{^}{L}}_{j}^{\left(k\right)}\left(\tau \right),$

where ${\stackrel{^}{L}}_{j}^{\left(k\right)}\left(\tau \right)=\prod _{l=1,l\ne j}^{{M}_{k}+1}\frac{\tau -{\tau }_{l}^{\left(k\right)}}{{\tau }_{j}^{\left(k\right)}-{\tau }_{l}^{\left(k\right)}}$, $1\le L\le {M}_{k}+1$, and ${M}_{k}$ is the number of the collocation points in interval ${S}_{k}$.

Integrating the expression in Eq. (18), Eq. (20) is expressed as:

20
${\stackrel{^}{X}}_{i+1}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{j}^{\left(k\right)}\right)={X}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{j}^{\left(k\right)}\right)+\frac{{t}_{f}-{t}_{0}}{2}\sum _{l=1}^{{M}_{k}}{I}_{jl}^{\left(k\right)}f\left({X}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{j}^{\left(k\right)}\right),{Z}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{j}^{\left(k\right)}\right),t\left({\stackrel{^}{\tau }}_{j}^{\left(k\right)},{t}_{0},{t}_{f}\right)\right),$

where ${I}_{ji}^{\left(k\right)}\left(j,l=1,2,\cdots ,{M}_{k};k=1,2,\cdots ,N\right)$ is a ${M}_{k}×{M}_{k}$ LGR integral matrix of LGR collocation point in $\left({\stackrel{^}{\tau }}_{1}^{\left(k\right)},{\stackrel{^}{\tau }}_{2}^{\left(k\right)}\cdots ,{\stackrel{^}{\tau }}_{{M}_{k}}^{\left(k\right)}\right)$:

21
${I}^{\left(k\right)}={\left[{D}_{2}^{\left(k\right)}\cdots {D}_{{M}_{k}+1}^{\left(k\right)}\right]}^{-1},{I}^{\left(k\right)}{D}_{1}^{\left(k\right)}=1.$

So, the relative and absolute error can be defined according to $X\left({\stackrel{^}{\tau }}_{l}^{\left(k\right)}\right)$ and $\stackrel{^}{X}\left({\stackrel{^}{\tau }}_{l}^{\left(k\right)}\right)$. The relative error ${E}_{i}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{l}^{\left(k\right)}\right)$ and absolute error ${e}_{i}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{l}^{\left(k\right)}\right)$ of the state variables in interval ${S}_{k}$ can be expressed as:

22
$\left\{\begin{array}{l}{E}_{i}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{l}^{\left(k\right)}\right)=\left|{\stackrel{^}{X}}_{i}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{l}^{\left(k\right)}\right)-{X}_{i}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{l}^{\left(k\right)}\right)\right|,\\ {e}_{i}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{l}^{\left(k\right)}\right)=\frac{{E}_{i}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{l}^{\left(k\right)}\right)}{1+\underset{j\in \left[1,2,\cdots ,{M}_{k}+1\right]}{max}\left|{X}_{i}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{j}^{\left(k\right)}\right)\right|},\end{array}\right\$

where $l=1,2,\cdots ,{M}_{k}+1$,$i=1,2,\cdots ,{n}_{x}$. And ${n}_{x}$ is number of the state variable.

Then the maximum relative error can be defined as:

23
${e}_{max}^{\left(k\right)}=\underset{\begin{array}{l}l\in \left[1,2,\cdots ,{M}_{k}+1\right]\\ i\in \left[1,2,\cdots ,{n}_{x}\right]\end{array}}{max}{e}_{i}^{\left(k\right)}\left({\stackrel{^}{\tau }}_{l}^{\left(k\right)}\right).$

### 4.2. Error range determination

Relationship between the estimated value $\left(\stackrel{^}{y},\stackrel{^}{u}\right)$ of LGR discrete interval and its analytical solution $\left({y}^{*},{u}^{*}\right)$ is:

24
${‖\left(\stackrel{^}{y}-{y}^{*}\right)‖}_{\infty }+{‖\stackrel{^}{u}-{u}^{*}‖}_{\infty }\le \frac{{c}_{1}{h}^{q}}{{N}^{q-\frac{5}{2}}},$

where ${c}_{1}$ is a constant; $N$ is number of the collocation point; $h$ is the interval width; $q$ is the minimum value of the order of consecutive derivatives solved and $N$; ${‖•‖}_{\infty }$ is interval norm.

### 4.3. Mesh refinement method

The solution result of Section 4.1 is used to determine whether the maximum relative error exceeds the tolerance error $\epsilon$ or not. If there is a grid relative error exceeding the tolerance error$\epsilon$, then the grid will be divided into smaller grids or the number of collocation point can be increased within the grid. When the state and control variables in the mesh are less smooth, the solution accuracy can be improved by refining the mesh. At the same time, when the state and control variables are smooth enough, the number of collocation points in the mesh is increased. Then the solution accuracy is improved by increasing the order of the interpolation polynomial.

### 4.3.1. Unsmooth mesh determination

It is assumed that the optimal solution on the $m$th interval, i.e. ${S}_{m}$, has been solved using the Radau pseudo-spectral method. And it is set that the solution of the state variable in $M$th iteration is ${X}^{\left(m\right)}\left(\tau \right)=\left[{X}_{1}^{\left(m\right)}\left(\tau \right),\cdots ,{X}_{{n}_{y}}^{\left(m\right)}\left(\tau \right)\right]$. Then continuous trajectory of the state and the control variable can be obtained by operating Lagrange interpolation on the state variable in ${S}_{m}$. So, the curvature of the trajectory can be expressed as:

25
${k}_{i}^{m}\left(\tau \right)=\frac{\left|{x}_{{m}_{i}}^{\text{'}\text{'}}\left(\tau \right)\right|}{{\left|1+{x}_{{m}_{i}}^{\text{'}}\left(\tau {\right)}^{2}\right|}^{\frac{3}{2}}},m=1,2,\cdots ,{M}_{k},i=1,2,\cdots ,{n}_{x},$

where ${k}_{i}^{m}\left(\tau \right)$ is the curvature of the $i$th state variable in ${S}_{m}$. Since each state variable has a different curvature change in ${S}_{m}$, the maximum value is taken as a judgment to determine whether the trajectory is smooth or not:

26
${k}^{m}\left(\tau \right)=\underset{i=1,2,\cdots ,{n}_{x}}{max}\left\{{k}_{i}^{m}\left(\tau \right)\right\}.$

If the curvature is satisfied with:

27
${k}^{m}\left(\tau \right)\ge \stackrel{-}{k},$

where $\stackrel{-}{k}$ is the set curvature value, indicating that the trajectory is relatively smooth in the interval, and the interval needs to be further divided. Otherwise, the trajectory is smooth in the interval, and the accuracy of the solution can be improved by increasing the number of points.

For control variables, the same method can be used to define the changes of the control variable.

### 4.3.2.1. Determination of sub-grids number

It is assumed that the curvature in grid ${S}_{m}$ in the $m$th iteration meets Eq. (27), and it is known from the above analysis that the mesh needs to be refined. Taking equal sign for Eq. (24), then Eq. (28) can be obtained:

28
${e}_{m}^{\left(M\right)}=\frac{{c}_{1}{\left[{h}_{m}^{\left(M\right)}\right]}^{q}}{{\left[{N}_{m}^{\left(M\right)}\right]}^{q-5/2}}.$

In order to obtain the desired minimum relative error $\epsilon$ in the $\left(M+1\right)$th iteration, Eq. (29) is obtained:

29
$\epsilon =\frac{{c}_{1}{\left[{h}_{m}^{\left(M+1\right)}\right]}^{q}}{{\left[{N}_{m}^{\left(M+1\right)}\right]}^{q-5/2}}.$

Since the grid ${S}_{m}$ has been refined, let ${N}_{m}^{\left(M+1\right)}={N}_{m}^{\left(M\right)}$, that is, the interval of ${S}_{m}^{M}$ in the $m$th interval and the iteration ${S}_{m}^{M+1}$ in $\left(M+1\right)$th iteration have the same number of collocation points. Eq. (30) can be obtained by combining Eq. (28) and Eq. (29):

30
$H=\frac{{h}_{m}^{\left(M+1\right)}}{{h}_{m}^{\left(M\right)}}={\left(\frac{{e}_{k}^{\left(M+1\right)}}{\epsilon }\right)}^{1/q},$

where $H$ is the number of sub-grids to be divided at the next step. The relative error in $\left(M-1\right)$th iteration is:

31
${\epsilon }_{k}^{\left(M-1\right)}=\frac{{c}_{1}{\left[{h}_{m}^{\left(M-1\right)}\right]}^{q}}{{\left[{N}_{m}^{\left(M-1\right)}\right]}^{q-5/2}},$

$q$ can be solved by combining Eq. (28). And then the number of sub-intervals can be determined according to Eq. (30). Eq. (24) demonstrates that the number of sub-intervals newly divided is at least $H$ to ensure that the relative error in the iterative process can meet the requirements. But for a specific problem, $H$ may be larger, so this will lead to a rapid increase in the amount of calculations, so it is necessary to set upper limit of $H$. In the study of the paper, the number of the largest subinterval of ${S}_{m}$ is obtained based on the ratio of relative and set error:

32
${H}_{max}=\left[\mathrm{l}\mathrm{o}{\mathrm{g}}_{N}\left(\frac{{e}^{\left(m\right)}}{e}\right)\right],$

where $\left[•\right]$ indicates taking an integer down. Eq. (32) gives the range of the largest subinterval. Its division process is as follows. When ${e}^{\left(m\right)}\to \epsilon$, Eq. (32) may cause the value of ${H}_{max}$ to be too large, resulting in the sub-grid being redundancy. When ${e}^{\left(m\right)}\to \epsilon$, ${H}_{max}$ decreases until it is 0 gradually. So ${H}_{max}$ always provides a reasonable number of subintervals. Then ${S}_{m}$ will be divided into:

33
$S=min\left(\left|H\right|,{H}_{max}\right).$

Eq. (33) gives the number of divided meshes but does not give the meshing position.

### 4.3.2.2. Determination of sub-grid position

In order to determine the location of the sub-grid, density function and cumulative distribution function are defined firstly. Let $f:\left[x,y\right]\to {R}^{+}\left(x,y\in R\right)$ be a non-negative Lebesgue integrable function, and ${\int }_{x}^{y}f\left(t\right)dt=1$, as well as $f$ be a density function. Then the integral of the density function on the defined domain is defined as the cumulative distribution function:

34
$F\left(t\right)={\int }_{x}^{y}f\left(\tau \right)d\tau .$

If $L$ collocation points ${\left\{{\tau }_{i}\right\}}_{i=1}^{L}=1$(${\tau }_{1}={T}_{k-1},{\tau }_{L}={T}_{k}$) are inserted in interval $\left[{T}_{k-1},{T}_{k}\right]$, and if position of ${\tau }_{i}$ is determined, then the position of ${\tau }_{i+1}$ is:

35
$F\left({\tau }_{i+1}\right)-F\left({\tau }_{i}\right)=\frac{1}{L}.$

The density function based on the curvature has a dense distribution at the position where the curvature is large, and the node distribution is sparse in a small curvature. Therefore, the cumulative distribution function can be used to determine the position of the newly added grid point. The density function can be defined by Eq. (25):

36
$\rho \left(\tau \right)={c}_{2}k{\left(\tau \right)}^{1/3},$

where ${c}_{2}$ is a selected constant and meeting the following constraint:

37
${\int }_{-1}^{1}\rho \left(\tau \right)d\tau ={\int }_{-1}^{1}{c}_{2}k{\left(\tau \right)}^{1/3}d\tau =1.$

Then the cumulative distribution function can be obtained:

38
$F\left(\tau \right)={\int }_{-1}^{\tau }\rho \left(\zeta \right)d\zeta ,\tau \in \left[-1,1\right].$

The position of the sub-grid is determined according to Eq. (33) and Eq. (34):

39
$F\left({\tau }_{i}\right)=\frac{i-1}{S},1\le i\le S+1.$

It is assuming that the solution in the grid ${S}_{m}$ does not satisfy the set error and Eq. (27). The default trajectory is smooth on interval ${S}_{m}$. Collocation points should be added to increase the order of the polynomial and to reduce its solution error. It is set that ${e}_{m}^{\left(M\right)}$ denote an error of the $M$th iteration on the interval ${S}_{m}$, and the grid width in the $\left(M+1\right)$th iteration is the same as in the $M$th iteration, i.e. ${h}_{k}^{M}={h}_{k}^{\left(M+1\right)}$. Then ${N}_{m}^{\left(M+1\right)}$ can be solved as:

40
${N}_{m}^{\left(M+1\right)}={N}_{m}^{\left(M\right)}{\left(\frac{{e}_{m}^{\left(M\right)}}{\epsilon }\right)}^{1/\left(q-5/2\right)},$

where ${N}_{m}^{\left(M\right)}$ is the number of collocation points of the $M$th iteration on grid ${S}_{m}$. Since the number of collocation points is an integer, Eq. (40) should be transformed into:

41
${N}_{m}^{\left(M+1\right)}=\left[{N}_{m}^{\left(M\right)}{\left(\frac{{e}_{m}^{\left(M\right)}}{\epsilon }\right)}^{1/\left(q-5/2\right)}\right].$

In order to ensure that the order of the polynomial on grid ${S}_{m}$ is not too high, ${N}_{max}$ is used to limit the order of the polynomial. If ${N}_{m}^{\left(M+1\right)}>{N}_{max}$, the number of collocation points on the grid is ${N}_{max}$.

### 4.4. Mesh refinement and collocation process updating

The process of grid refinement and collocation point number optimization is as follows:

Step 1. Setting the initial grid and the number of collocation points in each grid.

Step 2. Resolving the nonlinear planning problem shown in Eq. (18) using the SNOPT software package.

Step 3. If the relative error resolved is less than the set relative error $\epsilon$, the calculation is terminated, and the result is outputted; otherwise, Step 4 is to be done.

Step 4. Determining whether the curvature of the trajectory meets the set value $k$, if it does not meet the requirement, go to step 4.1; otherwise, go to step 4.2.

Step 4.1. Dividing the grids according to the method in Section 4.3.2 and setting the number of collocation points.

Step 4.2. Adding the number of collocation points according to the method in Section 4.3.3 or dividing the grid.

Step 5. After all the grids and collocation points have been updated, the calculation result of previous step is taken as the initial value, and then go to Step 2.

## 5. Numerical simulations and experimental verification

### 5.1. Numerical simulations

In order to verify the effectiveness of the method proposed in this paper, the simulations are made in MATLAB. They show that the vehicle parameters are according to Ref. .

### 5.1.1. Double Lane Changing Condition

In the test of driving stability of ordinary vehicles as well as in the test of tracking ability of intelligent vehicles, the double lane changing condition is a test section used frequently. However, since the track function of standard double lane changing condition is a piecewise function with poor continuity, in order to restore the true driving track under the double lane changing condition, the hyperbolic function is used to perform nonlinear fitting. And the following results are obtained:

42
$Y\left(X\right)=\frac{4.05}{2}\left(1+\mathrm{t}\mathrm{a}\mathrm{n}\left(a\right)\right)-\frac{5.7}{2}\left(1+\mathrm{t}\mathrm{a}\mathrm{n}\left(b\right)\right),$
43
$a=\frac{2.4}{25}\left(X-27.19\right)-1.2,$
44
$b=\frac{2.4}{21.95}\left(X-56.46\right)-1.2,$

where $X$ and $Y$ are the longitudinal and lateral positions respectively.

Fig. 2 shows the results of the lateral displacement, steering wheel angle, lateral acceleration, yaw rate and longitudinal speed under the double lane change road condition with the initial speed of 108 km/h.

From Fig. 2(a), it can be seen that the vehicle can track the double lane change road basically with the control of the proposed algorithm.

From Fig. 2(b), it can be seen that a large steering wheel angle is required to enter and exit the corners of the double lane change road. Especially, when entering the first corner and exiting the fourth corner of the double lane change road, the steering wheel angle reaches its maximum values.

From Fig. 2(c), it also can be seen that the lateral acceleration is higher at 60 m, 100 m, 430 m and 490 m; that indicates that when entering and exiting the corners of the double lane change road, the vehicle will generate the highest lateral acceleration.

From Fig. 2(d), it also can be seen that the yaw rate is higher at 60 m, 100 m, 430 m and 490 m; that indicates that when entering and exiting the corners of the double lane change road, the vehicle will generate a greater yaw rate.

From Fig. 2(e), it also can be seen that to complete the process of passing the double lane change road with minimum time, the vehicle must establish continuous acceleration.

Fig. 2Simulation results of the state variables under the double lane change condition a) Lateral displacement b) Steering wheel angle c) Lateral acceleration d) Yaw rate e) Longitudinal speed

### 5.1.2. Overtaking condition

Overtaking is a common operation mode in the process of vehicle driving.

Figs. 3(a)-(d) show the simulation results under the overtaking condition.

From Fig. 3(a), it can be seen that the vehicle can complete the processing of overtaking basically with the control of the proposed algorithm.

From Fig. 3(b), it can be seen that a large steering wheel angle is required to enter the adjacent lane from the current lane. Especially, in order to maintain stability after entering the adjacent lane, a large steering wheel angle is still required.

From Fig. 3(c), it also can be seen that a high lateral acceleration is required when entering the adjacent lane from the current lane.

From Fig. 3(d), it also can be seen that the yaw rate is greater at 15 m, 36 m; that indicates that when entering and exiting the current lane, the vehicle will generate a greater yaw rate.

Fig. 3Simulation results of state variables under overtaking condition a) Lateral displacement b) Steering wheel angle c) Lateral acceleration d) Yaw rate

### 5.2. Accuracy and efficiency verification

A comparison is made in order to verify the superiority of this method over the AMRC and GPM methods.

### 5.2.1. Double lane changing condition

Figs. 4(a)-(c) show the simulation results of the steering wheel angle and the error of steering wheel angle obtained by the AMRC and GPM methods in the double lane changing condition.

From Fig. 4(a), it can be seen that the trend of the curves of the steering wheel angle obtained by AMRC and GPM methods is basically consistent. However, it can be seen from Fig. 4(b) that the error of steering wheel angle calculated by the AMRC is smaller than that of the GPM method. This is because the method in this paper has fewer collocation and discrete equations and sparse Jacobian matrix, so it improves the speed of solving NLP and reduces the time of optimization. Fig. 4 indicates that under the same condition, it is required a less number of collocations and meshes, so the calculation accuracy is higher. This ensures the superiority of the proposed method. And also, it can be seen from Fig. 4(c) that the mean value of estimation error of steering wheel angle of the GPM method is 0.06 deg, while the one of the AMRC method is 0.028 deg, which is smaller than that of the GPM method. This ensures the higher calculation accuracy of the proposed method.

### 5.2.2. Overtaking condition

In order to analyze the computational efficiency of the proposed method, the calculation accuracy is compared between the proposed method and the GPM method. The comparison results are shown in Table 3. It can be seen from the table that for the same calculation accuracy, the number of iterations and the calculation time of the hybrid method were lower than that in the direct multiple shooting method. This demonstrated that the proposed hybrid method had better solving efficiency.

Fig. 4Comparison results between AMRC and GPM under double lane changing condition a) Steering wheel angle b) Error of steering wheel angle c) Mean value of estimation error of steering wheel angle

Table 3Comparison of number of iterations and calculation time of each method.

 Calculation accuracy (m) Number of iterations Time of calculation (s) AMRC GPM AMRC GPM 10-2 2 4 27.36 69.47

### 5.3. Experimental verification

A virtual test adopting the CarSim software is conducted to verify the feasibility of the simulated results.

CarSim vehicle dynamics software can quickly build a vehicle dynamics model and driving scene of the simulated vehicle by configuring parameters, and can also conduct a real-time simulation. The vehicle dynamics model is developed in the CarSim, which is a professional vehicle system simulation software. It is adopted by many automobile manufacturers and parts suppliers in the world, and it is the standard software of the vehicle industry. CarSim adopts a system-oriented parametric modeling method. Users only need to select each vehicle component module from the model database and configure the corresponding parameters as required to complete quickly the construction of the vehicle dynamics model. CarSim offers a variety of vehicle models including sedans, pickups, and SUVs. Carsim is the vehicle state research software that uses the Vehicle Sim technology based on the vehicle dynamics and takes the simulation as the technical method. During the project development and research, the use of Carsim can greatly reduce the time spent on development, testing and planning, and allow proposing better vehicle dynamics decision-making methods. For the research projects, such simulation technology can theoretically prove the expected results, which can provide a theoretical basis for a later research and development.

In this paper, the vehicle dynamics model of the driving simulation system is quickly built in the CarSim by configuring the parameters, and the control platform and motion simulation platform are developed by using its data interface. This improves the complexity of vehicle dynamics modeling in the process of the driving simulation system development, can reduce the difficulty of vehicle dynamics modeling in the driving simulation system, and simplify the development process.

From Fig. 5, it can be seen that for the lateral displacement and the steering wheel angle as well as longitudinal speed, there are discrepancies between the simulated and the virtual test values. This is because that the model in the virtual test ignores the nonlinearity of the steering system and suspension system. However, the trend of the curves in Figs. 5(a)-(c) is enough to illustrate the correctness of the proposed method.

Fig. 5Experimental results of lateral displacement and steering wheel angle as well as longitudinal speed a) Lateral displacement b) Steering wheel angle c) Longitudinal speed

## 6. Conclusions

To solve the problem of minimum time lane changing problem of vehicle handling inverse dynamics, this paper proposes an adaptive mesh refinement and collocation optimization method on the basis of traditional pseudo-spectral method and based on a nonlinear 7-DOF vehicle dynamics model as well as the Magic formula tire model. The simulation results show that the adaptive mesh refinement and collocation optimization method proposed in this paper can solve the problem of minimum time lane changing problem of vehicle handling inverse dynamics with high accuracy by approximating the optimal solution by densifying the mesh and adding collocation nearby. In the process of optimization solution, the method proposed in this paper can adaptively increase the number of collocations and reduce the number of meshes which can generate fewer total collocation in the solution area reducing the number of nonlinear programming equations after discretization, and improving the solving speed of the SQP method.

28 November 2022
Accepted
09 May 2023
Published
01 August 2023
SUBJECTS
Vibration in transportation engineering
Keywords
optimal control
vehicle lane changing
minimum time
mesh refinement
collocation method
Acknowledgements

This research was supported by the Science and Technology Program Foundation of Weifang under Grant 2015GX007. And also, this research was financially supported by the Open Research Fund from the State Key Laboratory of Rolling and Automation, Northeastern University, under Grant 2021RALKFKT008. The first author gratefully acknowledges the support agency.

Data Availability

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Author Contributions

Yingjie Liu: mathematical model and the simulation techniques. Dawei Cui: spelling and grammar checking as well as virtual validation. Wen Peng: virtual validation.

Conflict of interest

The authors declare that they have no conflict of interest.