Two-dimensional vehicular movement modelling at intersections based on optimal control

Modeling traffic flow at intersections is essential for the design, control, and management of intersections. A challenging feature of microscopic modeling vehicular movement at intersections is that drivers can choose among an infinite number of alternative traveling paths and speeds. This makes it fundamentally different from structured straight road sections with lanes. This study proposes a novel method to model the trajectories of vehicles in two-dimensional space and speed. Based on optimal control theory, it assumes drivers schedule their driving behavior, including the steering and acceleration, to minimize the predicted costs. The costs contain the running costs, which consist of the travel time and driving smoothness (longitudinally and laterally), and the terminal cost, which penalizes the deviations from the desired final state. Different than conventional methods, the vehicle motion dynamics are formulated in distance rather than in time. The model is solved by an iterative numerical solution algorithm based on the Minimum Principle of Pontryagin. The descriptive power, plausibility, and accuracy of the proposed model are investigated by comparing the calculated results under several cases, which can be solved from symmetry or analytically. The proposed model is further calibrated and validated using empirical trajectory data, and the quality of the predicted trajectory is confirmed. Qualitatively, the optimal trajectory changes in the range of the shortest path and smoothest path under different weights of the running cost. The proposed model can be used as a starting point and extended with more considerations of intersection operation in the real world for future applications.

To describe the vehicular movement in the inner area of the intersection, the description model should be twodimensional. Several works have been done using various methods based on the granular flow concept (Helbing et al., 2005;Hoogendoorn and Daamen, 2005). One approach is to extend the traditional one-dimensional car-following model to a two-dimensional model. Nakayama et al. (2005) constructed a two-dimensional optimal velocity model for pedestrians. Pedestrians are treated as identical particles moving in the two-dimensional space, and each particle decides its optimal velocity depending on distances to other particles. On this basis, Xie et al. (2009) proposed a twodimensional car-following model to describe the mixed traffic flow at unsignalized intersections. The interaction between a vehicle and other motorized/non-motorized vehicles were considered. Another approach is the social forcebased model, which can be applied to describe the pedestrian flow (Zeng et al., 2014), cyclist flow (Huang et al., 2017) and vehicular flow (Ma et al., 2017b) at intersections. Ma et al. (2017a) proposed a three-layered "plan-decisionaction" model to simulate the moving of turning vehicles at mixed-flow intersections and further used the social force model to generate vehicle movements in the operation layer (Ma et al., 2017b). Moreover, several social force-based vehicle moving models were established to simulate the vehicular traffic flow with the consideration of no lane division (Fellendorf et al., 2012;Huynh et al., 2013;Yang et al., 2018). The trajectory is obtained by the comprehensive result of several types of forces, such as the self-driven force, the repulsive force, the attractive force, and other forces for specific conditions. Moreover, the results of the optimal trajectory guidance for the autonomous vehicle also describe the vehicular movement at intersections. Numerous works considered the optimal trajectory guidance for autonomous vehicles traversing crossing the intersection. Some studies focused on establishing cooperative vehicle intersection control algorithms to enable autonomous vehicles to pass the intersection cooperatively without traffic signals (Ahmane et al., 2013;Lee and Park, 2012;Yu et al., 2019), while some others aimed to optimize the vehicular trajectory under the consideration of the signal control (Jiang et al., 2017;Kamalanathsharma et al., 2015) or to optimize of the trajectory and signal control in a unified framework Liu et al., 2019;Xu et al., 2018;Yu et al., 2018). Dresner and Stone (2004) propose a reservation-based system at intersections under the condition of fully autonomous vehicles. The system aims to maximize the efficiency of moving cars through intersections with 0 probability of collisions. The reservation-based method was further enhanced to improve various operational issues (Dresner and Stone, 2008;Li et al., 2013). Medina et al. (2015) developed a cooperative intersection control method by defining a virtual intervehicle distance between vehicles driving on different lanes. Based on the virtual platooning concept, the model can guide the vehicles to cross the intersection safely. Zohdy and Rakha (2016) developed a system, namely iCACC (intersection cooperative adaptive cruise control), to optimizes the movement of CACC equipped vehicles with the objective of minimizing the intersection delay while ensuring operational safety. Li et al. (2019) proposed a speed guidance algorithm for automated vehicles to ensure safety and operation efficiency at intersections. The trajectories of automated vehicles were formulated in a 3-dimensional way, in which the vehicle width and length and time continuity were considered. Guo et al. (2019) further proposed an algorithm to optimize the trajectories of automated vehicle and signal control simultaneously for best system performance in terms of travel time, fuel consumption, and safety. However, in these models, the vehicular trajectories are assumed to be along the predetermined traffic lanes. Bichiou and Rakha (2018) developed a novel algorithm is to obtain an optimal solution for each possible scenario a vehicle may encounter while crossing an intersection. The detailed motion of the vehicle and the constraints on avoiding collisions were considered (Fadhloun and Rakha, 2019;Rakha et al., 2001;Rakha et al., 2009). The promising effectiveness of the algorithm in reducing delay and emissions was analyzed by comparing it with other intersection control strategies. The algorithm was further simplified to make it practical for real-time implementations (Bichiou and Rakha, 2019). These models also use time as the main independent variable. The control inputs are the rate of change of the velocity and its vector orientation, which are appropriate for vehicle control. However, drivers cannot change the angular velocity of the vehicle directly, but they can adapt their steering wheel angle, which changes the curvature of the trajectory.
Compared to the urban street segments or freeways, a challenging feature of microscopic modeling vehicular movement at the inner area of intersections is that drivers can choose among an infinite number of alternative running paths and speeds according to their desires. The traditional microscopic models are lane-based (Munigety and Mathew, 2016;Rahman et al., 2013), i.e., vehicles should travel in a pre-determined traffic lane or move from one lane to another, which may not accord with the fact that the concept of the traffic lane is significantly weakened at the inner area of intersections. Some newly proposed methods, including the two-dimensional car-following model, the social force-based model, and the optimal trajectory guidance model, have the possibility of describing the two-dimensional moving process of vehicles. However, the trajectory is not obtained directly from the behavior of drivers, such as turning the steering wheel and pushing the brake or throttle pedals. Therefore, it is still a challenge to make a realistic description of the vehicular trajectory driven by the human driver at intersections. The motivation of the paper is to model the trajectory of the manual driving vehicle at intersections. A twodimensional vehicular movement model is proposed based on optimal control. In methodology, the study incorporates both steering and acceleration as model inputs and chooses the travelled distance as the main independent variable instead of a more intuitive variabletime. In this way, the longitudinal and lateral controls are decoupled. That means that a particular heading angle as a function of the distance leads, independent of the acceleration function, to a particular path in space. Therefore, the model resembles human driving inputs, consisting of steering (steering wheel) and longitudinal acceleration (pedals). The model is solved by an iterative numerical solution algorithm based on the Minimum Principle of Pontryagin. The descriptive power, plausibility, and accuracy of the proposed model are validated by comparing the calculated results for several cases with expected or analytically calculated results. Determining the exact trajectory can help the researchers and engineers to evaluate the intersection design and control scheme more accurately. For example, based on the proposed model, the distribution of the moving trajectory can be obtained under different design schemes. Furthermore, the proposed model can be used in the microscopic simulation software to improve the reality of the result, instead of predetermining the running path.
The article is organized as follows. Section 2 establishes the utility-based two-dimensional vehicular movement model. Section 3 analyzes two specific conditions, namely the fastest path and the smoothest path. Section 4 validates the descriptive power, plausibility, accuracy, and the quality of the predicted trajectory of the proposed model. Section 5 discusses the expansion of the model with some considerations of the operation in the real world. Section 6 summarizes the research results and discusses future research directions.

Optimal vehicular movement model formulation
At intersections, vehicles move from an approach lane to an exit lane, as shown in Fig. 1. In continuous space, an infinite number of paths from origin to destination is possible. In this research, we assume that the driver will minimize the predicted disutility/cost while simultaneously choosing the path and speed. The remainder of this section describes this problem; we first define the modelling variables, then the cost functions, followed by constraints and the solution algorithm.

Vehicle motion dynamics
We define state as a function of the distance traveled ( ) by the subject vehicle; the state contains the xcoordinate ( ), y-coordinate ( ), the heading angle ( ), and the pace ( ), i.e., the reciprocal of the vehicular velocity, as shown in Eq. (1). Choosing the travelled distance ( ) as the main independent variable compared to a more intuitive variabletime has considerable advantages. First, the longitudinal and lateral controls are decoupled. That means that a particular steering function as function of leads, independent of the acceleration function, to a particular path in space. Moreover, within this lateral movement, the concept of the is equal to the curvature of the path of the vehicle.
where, denotes the state of a vehicle at moving distance ; and are the plane coordinate of a vehicle at moving distance , in m; is the heading angle at moving distance , in rad; is pace, i.e., the reciprocal of the velocity of a vehicle at moving distance , in s/m.
In this study, it is assumed that the initial and desired terminal states are given. The initial state of a vehicle 0 is described with: (2) while the desired terminal state of the vehicle is described with Eq. (3), in which the terminal pace is free: To describe the trajectory choice behavior, the driver uses a utility model to estimate and predict the trajectory costs. Then, the driver will determine his future state using the state prediction model (the system dynamics), which can be expressed by Eq. (4). According to Eq. (4), and constitute the control , as shown in Eq. (5). The control variable can be seen as the control of the steering wheel, while another control variable can be seen as the control of the brake or throttle pedals. In this way, we control the vehicle as the human driver does.
where, is the curvature (reciprocal of the turning radius) of a trajectory at moving distance , in rad/m; is the parameter indicating the acceleration of a vehicle at moving distance , in s/m 2 , a positive value indicates decelerating, a negative value indicates accelerating.
where, denotes the control vector of a vehicle at moving distance .

Moving cost
Consider the control for a trajectory from the initial state 0 to the actual terminal state . Given the control function ( ), the state of the vehicle along the full trajectory can be determined, which will also yield the disutility (cost), as shown in Eq. (6).
In this equation, denotes the terminal cost, which reflects the costs at the end of the trajectory; denotes the moving distance at the end of the trajectory; denotes the actual terminal state; denotes the running cost, which reflects the costs occurring during a trajectory. We will now elaborate on these costs.

Composition of terminal cost
The aim of the terminal cost is to get a trajectory close to the desired position and heading angle. In this line, the terminal cost can be specified as Eq. (7) where deviations of the actual final state (indicated with subscript ) from the desired final state (indicated with subscript ), induce cost: In this equation, denotes the relative weights of the terminal cost, which can be set as a relatively large number to ensure the vehicles reach the desired final state (position and angle).

Composition of running cost
The running cost reflects the different cost aspects considered by the drivers during driving. We model it as a sum of different elements: = ∑ where, denotes the contribution of the cost aspect ; denotes the relative weights of the running cost aspect . We will consider 3 running cost terms. Firstly, the travel time of the trajectory is one of the contributing costs. This is formulated as Eq. (9). Since the pace is integrated (Eq. (6)) over the trajectory, this yields the travel time to the final position.

=
(9) The remaining two terms reflect the driving smoothness. Drivers might not prefer to drive with a sharp turn, nor prefer to drive with frequent acceleration and braking. The degree of the driving smoothness is defined as a quadratic function of the vehicular lateral acceleration (centripetal acceleration) and longitudinal acceleration, as shown in Eqs. (10) and (11), respectively. where, is the vehicular lateral acceleration, in m/s 2 ; is the vehicular longitudinal acceleration, in m/s 2 .

Constraints
To conform to reality, the running speed, the curvature, and the acceleration should be limited according to the traffic rule and the characteristic of a vehicle. The running speed should be restricted within a reasonable minimum and maximum speed limitation, and the reversing is prohibited, as shown in Eq. (12). The curvature of a moving path is bounded by a given minimum turning radius, as shown in Eq. (13). The other control variable that indicating the vehicular acceleration should also be in a reasonable range, as shown in Eq. (14). Please note that the proposed model aims at generating motions in 2-dimensional space and time with a minimal set of parameters that can be used to describe traffic operations at intersections. In regular human driving, the limits of the vehicular physical capabilities are not reached. Therefore, the simplified vehicle constraints on its motion are considered in our work; this might differ for automated driving.

Solution
In this section, we describe the solution method. An iterative numerical solution algorithm was used based on the Minimum principle of Pontryagin.

Minimum principle of Pontryagin
The utility optimization paradigm implies that the drivers will choose the control that can minimize the total cost, which is a function of the initial state, the desired terminal state, and the control. Let * denotes the optimal control to minimize the predicted cost, as shown in Eq. (15). Let * denotes the value function of the cost of a vehicle when the optimal control * is used, as shown in Eq. (16). * = arg min ( , ) (15) * = min ( , ) = ( , * ) The minimum principle of Pontryagin is used to solve the optimal control problem. The method entails defining the Hamiltonian ℋ, as shown in Eq. (17), which can be specified as Eq. (18).
where ℋ denotes the Hamiltonian; ′ is the transpose of , which denotes the so-called co-state or marginal cost of the state . These costs reflect the relative extra cost due to making a small change Δ on the state . 1 , 2 , 3 , and 4 are the co-states corresponding to the four elements of the state vector , , , and , respectively. Using the Hamiltonian, we can derive the following necessary condition for the optimal control * : ℋ( , , * , ) ≤ ℋ( , , , ) For the running cost function shown in Eq. (7), using the necessary condition (19) we can find the expression for the optimal control, as shown in Eqs. (20)-(21).
Moreover, for the co-state , we can determine the costate dynamics as shown in Eq. (22), subject to the terminal conditions at the end of the control horizon, as shown in Eq. (23). Then, we get specific co-state dynamics, as shown in Eq. (24), and the transversal conditions, as shown in Eq. (25).

Iterative numerical solution algorithm
To solve the model mentioned above, we present a numerical solution. This solution approach is based on iteratively solving the state dynamic equation forward in the distance, and subsequently solving the co-state equation backward in the distance Wang et al., 2012). The following procedure summarizes the algorithm: Step 1: Initialization. Input the initial and desired terminal states ( 0 and ), choose the adjustment factor (0.999 is set as the initial value, and it can be dynamically enlarged when the algorithm cannot converge), set the iteration number = 1, set the initial co-state 0 = 0, and set the control step.
Step 2: Solve the state dynamic equation forward in time according to Eqs.
Step 3: Solve the co-state dynamic equation backward in time according to Eqs. (22) and (23). Go to step 4.
Step 4: Update the effective co-state according to Eq. (26). Go to step 5.
Step 5: Stopping criteria. The iteration will be stopped if the difference of the co-states of the backward and forward is less than the threshold (0.1 is set), as shown in Eq. (27); otherwise, set = + 1, and go back to step 2.

Specific condition analysis
The iterative numerical solution algorithm is applicable for a wide range of conditions, but the result is an approximate one because the precision of the control step and convergence threshold is limited. In this section, we solve the proposed model analytically under some specific conditions. Two specific conditions will be analyzed, namely the fastest path and the smoothest path. The numerical algorithm results will be compared and verified with the analytical results in section 4. This will yield a face validation for the solution method.

The fastest path
In this case, the expected travel time is the only running cost considered. Therefore, the Hamiltonian function is specified as Eq. (28), and the co-state dynamics as Eq. (29). We can derive the necessary condition (19) for optimal control.
To minimize the Hamiltonian, for the curvature ( ), either one of the following two cases holds.
(1) Case one: curvature ℋ = 3 = 0. Then, . Therefore, the curvature should equal to 0 ( = 0). (2) Case two: ℋ ≠ 0. Then, in order to minimize the ℋ shown in Eq. (28), we should minimize 3 . Therefore, we should select the boundary of the curvature ( = ±  (25). This means 4 is always nonnegative ( 4 ≥ 0). The optimal control of should select the minimal value under the constraints Eqs. (12) and (14), as shown in Eq. (31). In summary, Eqs. (30)-(31) is the analytical solution of the proposed model under the condition of the fastest path. The conclusion is that the optimal control of the fastest path is the concatenation of the given maximum curvature ( = ± 1 min ) and fixed heading angle with the curvature equals to 0 ( = 0) with the maximum acceleration (the minimum value of ) under the restrictions of the speed and acceleration. This conclusion will be used to validate the plausibility of the numerical algorithm results in section 4.2 (see the red line in Fig. 3(g) and (i)). Moreover, based on the analytical solution, the accuracy of the numerical algorithm results will be verified in section 4.3 (see Table 3).

The smoothest path
In this case, the driving smoothness is the only running cost considered. Therefore, the Hamiltonian function turns to be Eq. (32), and the co-state dynamics turns to be Eq. (33). Then, using the necessary condition (19), the expression for the optimal control for the curvature ( ) and the acceleration ( ) can be drawn, as shown in Eqs. (34) and (35), respectively.
They are the same as Eqs. (20)- (21), which is very difficult to obtain a full analytical solution. However, we can still analyze some general results for the control variable of acceleration ( ). We get d 4 d ≥ 0 according to Eq. (33), because 2 , 3 , 2 , 2 , and are all non-negative. It means the value of 4 will not decrease. Meanwhile, we have the terminal value of 4 ( ) = 0 according to Eq. (25). Therefore, 4 should be no larger than 0 ( 4 ≤ 0 ).
Additionally, 3 and are larger than 0. Then, = − 4 6 3 is nonnegative, as shown in Eq. (36). Therefore, a partial analytical solution can be drawn that the optimal control of is nonnegative under the smoothest path condition, as shown in Eq. (36). The conclusion is that the optimal control of the smoothest path is decelerating or keeping the same speed. This conclusion will be used to validate the plausibility of the numerical algorithm results in section 4.2 (see the red line in Fig. 3

Model validation
The descriptive power, plausibility, accuracy, and the quality of the predicted trajectory of the proposed model will be validated in this section according to the following four steps: (1) validate the descriptive power of a solution to the model by enumerating common moving cases in the intersection; (2) validate the plausibility of the solution of model behavior by analyzing change tendency of the path, the two control variables ( and ), and the running costs under different weights of the running costs; (3) validate the accuracy of the solution algorithm by comparing the difference between the results of the proposed algorithm and the analytical solution; (4) calibrate the model parameters and validate the quality of the predicted trajectory using empirical data.

Descriptive power validation
To validate the proposed solution algorithm can be used to solve all conditions, the common moving cases in the intersection are enumerated in Table 1, which contains U-turn, left turns with various angles, through movements with/without an offset, and right turns with various angles. The other input parameters are shown in Table 2. The results are illustrated in Fig. 2, which is a manifestation that the proposed solution algorithm can generate optimal control results under all tested scenarios.

Plausibility of model behavior
For reasons of simplicity in explanation, the left-turn scenario (Scenario F in Table 1) will be used in the hereafter discussion. The path ( , ), the heading angle ( ), the pace ( ), the curvature control ( ), and the acceleration control ( ) are analyzed, as illustrated in Fig. 3. Two sub-scenarios are discussed in this section. Scenario 1 is a free terminal pace scenario, which is the same as Scenario F in Table 1. Scenario 2 is a restricted terminal pace scenario, in which the terminal pace is restricted to be equal to the initial one. The terminal state in Scenario 2 turns to be [10, 16, /2, 1/8]. The other input parameters are shown in Table 2.
Generally, the path (see Fig. 3(b)), the running state (see Figs. 3(d) and 3(f)), and the curvature and acceleration controls (see Figs. 3(h) and 3(j)) are symmetric under the restricted terminal pace scenario. It indicates that the trajectories are the same when driving from one point to another and backward when the initial state and the terminal state are the same, which is in accordance with the expectation. However, the curvature and acceleration controls, the running state, and the path are nonsymmetric when the pace is changeable. It is because of all the three aspects of the running cost ( 1 , 2 , and 3 ) are related to the pace. The two control variables, and , will adjust according to the pace state to obtain the minimum running costs.
For the running path, as we expected, it is always within the boundary of the fastest path and the smoothest path. The path becomes straighter and the turning becomes sharper with the increase of the weight of the travel time cost. Accordingly, the length of the running path decreases with the increase of the weight of the travel time cost. It is because when the weight of the travel time cost increases, the vehicle tends to select the fastest path (the black line); on the country, the vehicle tends to select the smoothest path (the red line) when the weight of the driving smooth cost increases.
For the change tendency of the heading angle, it changes more homogeneous with the increase of the relative weight of the unsmooth running cost. As the red line shown in Fig. 3 (c) and (d), it is almost straight, which indicates the heading angle changes uniformly. It is in accordance with the definition of the "smoothest path". On the contrary, the black line is the most unsmooth path. It is composed of two segments with the largest slope and one segment with a slope of 0. The two segments with the largest slope indicate the vehicle runs with the largest curvature, while the segment with a slope of 0 indicates the vehicle runs straight ahead. It can be further verified by Fig. 3 (g) and (h), where we can see the change of the curvature control along the path. The red line (smoothest path) is the most stable one. The black line (the fastest path) has the largest changes (maximum curvature or 0), which is in accordance with the analytical results of the fastest path shown in section 3.1.
For the change tendency of the pace (p), As expected, the changes of p under the fastest path is the most dramatic. The p of the fastest path (black line) decreases (acceleration in velocity) with a sharp slope then maintains at the minimum value (the maximum velocity limitation). It is in accordance with the change of the in Fig. 3 (i) and (j), where only the minimum value or the 0 of the is selected. It is in accordance with the analytical results of the fastest path shown in section 3.1. Moreover, from Fig. 3(j), we can find the value of always larger or equal to 0 for the smoothest path (red line), which is in accordance with the analytical results of the smoothest path shown in section 3.2. For other paths except for the fastest and smoothest paths, the change of the running pace depends on the weights of the travel time cost and unsmooth running cost. When the travel time cost makes a primary contribution to the total cost, the vehicle prefers to accelerate (negative value in , hence a reduction of pace and increase of speed)) to minimize the total cost. Fig. 4 shows the change in the travel time cost and the driving smoothness cost with the change of the weights. As expected, the travel time cost decreases and the driving smoothness cost increases with the increase of the weight of travel time cost.
In summary, the analysis shows that the proposed model is sensible. The vehicle dynamics are shown to have symmetry when the terminal pace is restricted to be equal to the initial one. Qualitatively, the optimal trajectory changes in the range of the shortest path and smoothest path under different weights of the running cost.

Accuracy validation
The result of the fastest path is further analyzed in this section by comparing the results drawn from the analytical solution in section 3 to validate the accuracy of the proposed algorithm. According to the input parameters shown in Table 3, the expressions of the real fastest path of the left turn (Scenario F in Table 1) are shown in Eq. (37). We can check the gap (the Euclidean distance for the path and the absolute error of the pace) between the calculated results based on the proposed algorithm and the real ones obtained from the analytical solution at each calculated step. The comparison results are shown in Fig. 5 and Table 3.
The results show that the differences are small. In this test case, the average error, maximum error, and RMSE (root-mean-square error) in path (position state) of the real fastest path are 0.04 m, 0.08 m, and 0.03 m, which very small compared to of the width of the traffic lane. Additionally, there is no observable difference in the pace. Therefore, the accuracy of the proposed algorithm is acceptable.

Validation with empirical data
In this section, the proposed model is validated using empirical data. The model parameters are firstly calibrated using the real vehicular trajectories collected in an intersection. Then, the calibrated model is validated using the data collected in another intersection.
Two intersections located in Shanghai, China, namely Youyi Road -Tieli Road and Zuchongzhi Road -Gaosi Road, are selected for parameter calibration and model validation, respectively, as shown in Fig. 6. The operation of the intersection is recorded by the unmanned aerial vehicle. The position of each vehicle at each frame (1/24 seconds) can be collected by using a specially developed video recognition software. Then, the real vehicular trajectory can be obtained.
(a) Intersection for parameter calibration (b) Intersection for model validation Fig. 6 Surveyed intersections The RMSE (root-mean-square error) of the trajectory errors at all time steps will be used as the indicator of the accuracy of the model. The trajectory error at each time step (1/24 seconds) is defined as the Euclidean distance between the generated trajectory and the real trajectory, as shown in Fig. 7(a).
For parameter calibration, the three relative weights of the running cost, 1 , 2 , and 3 , need to be calibrated in the model. Since they are relative weights, without loss of generality, 1 is set to be 1. Then, there are only two parameters to calibrate, = { 2 , 3 }. For each collected vehicle, the was calibrated using the enumeration method according to the following five steps. (1) Extract the initial and terminal states of the vehicle. (2) Enumerate all the from 0 to 0.5 with an interval of 0.001. (3) Generate the trajectory under each using the proposed model. (4) Calculate the corresponding fitness value (RMSE). (5) Obtain the optimal value of according to the fitness (smallest RMSE). Fig. 7(b) and (c) illustrate an example of the parameter calibration for left-turn and through vehicles, respectively. The results of the parameter calibration of 40 vehicles are shown in Table 4. One can find that the values of RMSE are smaller than 2.0 m for all the tested vehicles. The average RMSE of trajectory error is 0.923 m. It means the distance between the generated trajectory and the real trajectory at each time is approximately one meter. If only considering the path and neglecting the time and speed, the average RMSE of path error is 0.486 m. Given the initial and terminal states, the real trajectories can be well represented by adjusting the weights of the costs. It indicates that the proposed model has a strong performance in describing the moving trajectory of the manual driving vehicle.
For model validation, it was conducted by comparing the real vehicular trajectories and the predicted trajectories with the calibrated . For each vehicle, the accuracy of the model was tested according to the following four steps. (1) Extract the initial and terminal states of the vehicle. (2) Generate the trajectory under each calibrated in Table 4 using the proposed model. (3) Calculate the corresponding fitness value (RMSE). (4) Obtain the accuracy of the model and the optimal value of according to the fitness (smallest RMSE). Fig. 8(a) and (b) illustrate an example of the model validation for left-turn and through vehicles, respectively. The results of the model validation of 20 vehicles are shown in Table 5. The average RMSE of trajectory error is 1.237 m. Although it is a little larger than the parameter calibration cases, the accuracy is still acceptable. Please note this is the result of the trajectory using other driver's preferences. In application, numerous real vehicular trajectories can be collected to draw the distribution of the calibrated parameters ( ). Then, we can represent the distribution of the trajectories by running the model many times using different values of drawn from the distribution.

Model extensions
In this paper, we establish a utility-based two-dimensional model to describe the trajectory of the manual driving vehicle at intersections. The model is a simplified representation of real-world driving, of which realism can be extended with elements such as the uncertainty of the manual driving, the preference of the driving path, and the avoidance of stationary and moving obstacles. This section discusses these possibilities further. Since the model is to describe the trajectory of human driving, we prefer to reflect these phenomena as the result of the trade-off between different costs instead of setting constraints. In short: the framework of control actions (splitting longitudinal and lateral controls and applying a Pontryagin principle for optimization) is the same, but additional costs are included for these considerations.

Uncertainty of manual driving
For a human driver, even if we assume that the driver knows the optimal control, it is still difficult or even impossible to control a vehicle exactly what he wants. In the hierarchical framework of strategic, tactical, and operational driving, as defined by Hoogendoorn and Bovy (1999), the vehicle operation is not necessarily accurate. The actual control * contains the uncertainty of the manual driving. Moreover, the driver will realize this difference and re-plan its desired trajectory (in the tactical layer). To overcome this modelling shortcoming, we can add the noise to the moving dynamics to reflect the uncertainty of the manual driving, as shown in Eq. (38). Additionally, the idea of re-planning the optimal control at each discretized control period can be applied. It means a driver will plan the entire control strategy from the current state to the terminal state, but only the first period of the control decision will be applied. After that, the driver will re-plan the control strategy from the updated state to the terminal state. The discretized control period can be pre-determined with the unit of distance or time. In this way, the uncertainty of manual driving can be reproduced. The actual running trajectories change for different drivers, even their weights to the moving costs are the same. We model it as: where, * and * are the optimal control of and ; 1 and 2 represent the uncertainty of the human driver control and the resulting effects on the subject's kinematics.
Using the same case study as in section 4.2 (Scenario F in Table 1), the random variables 1 and 2 obey the normal distribution with the expectation of 1 = 2 = 0, and the standard deviation of 1 = 2 = 1. The discretized control period is set to be 4 m. The weights of the running cost 1 , 2 , and 3 are set to be 1, 0.001, and 0.01, respectively. Other parameters remain the same as those in Table 2. Fig. 9 illustrates the optimal control and the manual driving of each control period. One can find that the trajectories of the manual driving are not the same as the optimal control. There are some errors. However, at the beginning of each control period, the optimal control path will be re-planned. Then, the human driver adjusts the control according to the re-planed route and reaches the terminal. Please note the uncertainty of manual driving changes for different drivers. Therefore, the actual running trajectories change for different drivers, even the distribution of the random variables 1 and 2 are the same. As an example, Fig. 10 illustrates 5 manual driving trajectories with the same distribution of 1 and 2 .

Preference of the driving trajectory
Human drivers also may prefer to drive in some space instead of others according to the layout of the intersection, which can be called as the preference of the driving trajectory. It means the running cost in different areas is different. This might be lanes, or closer to the side. An additional running cost term reflecting this phenomenon can be added in this framework, which is related to the position (coordinates) of the vehicle, as shown in Eq. (39).
where, ( , ) is the cost function related to the coordinates.
The explicit formulation of Eq. (39) is related to the layout of the intersection and the specific geometric design element that planned to explore. As an example, we show a case with a turning movement, like introduced in section 4.2. In this case, we consider two sets of line-markings that are used to guide vehicles. One makes a smooth turn, and one a sharp turn (Fig. 11). The equations are as follows.
Assuming the curve of the guide line is given, the cost function can be defined as the quadratic function of the distance between the vehicular position ( , ) and the guide line curve ( ′, ′), as shown in Eq. (40).
( ′, ′) = 0 (41) As shown in Fig. 11, the two guide line curves can be expressed by Eqs. (43) and (44) for the smooth guide line and the sharp guide line, respectively.
We can derive the necessary condition (19) for the optimal control for both guide line scenarios. The results of the optimal control for the smooth and sharp guide line scenarios under various weights of the driving preference cost are shown in Fig. 12 and Fig. 13, respectively. The weights of the running cost 1 , 2 , and 3 are set to be 10, 0.001, and 0.01, respectively. The weights of the running cost 4 changes from 0 to 5. Other parameters remain the same as those in Table 2.
One can find that in both guide line scenarios, with the increase of the weight of the driving preference cost ( 4 ), the running trajectory becomes closer and closer to the designed guide line, while the value of the velocity almost remains the same. Comparing the two scenarios of the guide line, the sharp guide line is more difficult to approach. Drivers will not follow the sharp guide line unless a huge weight is given for the running cost of driving preference, which can be a traffic rule in reality. Moreover, we noticed that the running cost Eq. (48) is a piecewise continuous function. It leads to the Hamiltonian function (49) becoming a piecewise continuous function. The results show that the proposed iterative numerical solution algorithm can handle this condition.

Interference with stationary and moving obstacles
In the real world, drivers should avoid objects during driving. This can also be captured by adding a new running cost term. The generalized obstacles can be the separation facilities, the piers of the elevated highways and pedestrian overpasses, and other vehicles. They can be treated as stationary or moving obstacles. Since the drivers aim to stay clear of the obstacles, the running cost caused by the obstacles can be defined as a monotonically decreasing function of the distance between a vehicle and the obstacle, as shown in Eq. (51). Then, the Hamiltonian function can be specified as Eq. (52), and the co-state dynamics turns to be Eq. (53). The difference between the stationary and moving obstacles lies in that the position of the moving obstacle ( , ) will change along the time while the stationary obstacle position will not. Please note, since the state of a vehicle is described by using the moving distance as the dependent variable in the proposed model, the moving distance should be transformed to the time when updating the position of the moving obstacles, as shown in Eq. (54). In a quite short distance, Eq. (54) turns to Eq. (55) by using the Euler forward difference.
The difference between the stationary and moving obstacles lies in that the position of the moving obstacle ( , ) will change along the time while the stationary obstacle position will not. Please note, since the state of a vehicle is described by using the moving distance as the dependent variable in the proposed model, the time should be obtained from distance when updating the position of the moving obstacles, as shown in Eq. (54). In a quite short distance, Eq. (54) turns to Eq. (55) by using the Euler forward difference.
= ∫ 2 1 (54) Using the same case study as in section 4.2, two scenarios of obstacles are considered: one is a stationary obstacle located at the coordinate (4.0, 8.5); the other is a moving obstacle driving straight forward from the position (12, 10) to the west with a constant velocity of 4 m/s. The weights of the running cost 1 , 2 , 3 , and 5 are set to be 10, 10 -3 , 10 -2 , and 10, respectively. The parameter that reflecting the influence range of the obstacle ( ) is set to be 0.75. Other parameters remain the same as those in Table 2. The results of the optimal control for stationary and moving obstacle scenarios under various weights of the obstacle cost are shown in Fig. 14 and Fig. 15, respectively.
For stationary obstacles, vehicles should adapt their path in space to go around the obstacle and have a longer path, as shown in Fig. 14. With the increase of the parameter from 0 (without the consideration of the obstacle) to 0.5, 0.75, and 1, which reflects the influence range of the obstacle, the minimum distance between the running trajectory and the stationary obstacle changes from 0.84 m to 1.53 m, 2.17 m, and 2.84 m, respectively. In reality, the value of the parameter should be adjusted according to the width of the vehicle and the size of the obstacle. For moving obstacles, the framework allows two ways to avoid the moving obstacle: one is slowing down to yield the moving obstacle, the other is speeding up to pass the conflict point before the moving obstacle arriving. The choice of the optimal control strategy depends on the weights of different aspects of the running costs ( ), as shown in Fig.  15. The color of the trajectory indicates the running time. There is a conflict between the controlled vehicle and the moving obstacle if the running cost of the obstacles is not considered (Trajectory 1: 1 = 10, 2 = 10 −3 , 3 = 10 −2 , and 5 = 0). The conflict can be successfully avoided after considering the obstacles in the running cost. With the decrease of the weights of the smooth-running costs, the trajectory changes from a longer detour path with lower travel speed (Trajectory 2: 1 = 10, 2 = 10 −3 , 3 = 10 −2 , and 5 = 10) to a shorter path with sharp turning and higher travel speed (Trajectory 3: 1 = 10, 2 = 10 −4 , 3 = 10 −3 , and 5 = 10). Reducing the weights of the smoothrunning costs results in aggressive behavior.

Conclusions
This study proposed a novel vehicular movement model at intersections, based on the assumption that the drivers are utility maximizers. The drivers schedule their activities, including the steering and acceleration, to minimize the predicted costs. The problem is described by an optimal control framework. The terminal costs of the deviations from the desired final state are considered, and the running costs of the travel time and driving smoothness are considered, which reflects a trade-off between the utilities of time and comfort. The intra-driver variation mechanism can be reflected by the weights of the costs. The proposed model is solved by an iterative numerical solution algorithm based on the Minimum Principle of Pontryagin. The descriptive power, plausibility, and accuracy of the proposed model are validated by comparing the calculated optimal control results with the common senses and the analytical solution under several cases. To operationalize the proposed model, it is extended by adding more considerations of the prevailing traffic conditions, including the uncertainty of the manual driving, the preference of the driving trajectory, and the avoidance of stationary and moving obstacles.
The main contribution of the study is the description of vehicle dynamics at intersections where the drivers can choose the optimal control among an infinite number of alternative running paths and speeds. The two control variables used in the model, and , are chosen to represent how the drivers schedule their driving behavior, including turning the steering wheel and pushing the brake or throttle pedals. Choosing the distance along the trajectoryrather than timeas an independent variable moreover allows to decouple speed and path, which is useful for trajectory optimization. Qualitatively, the optimal trajectory changes in the range of the shortest path and smoothest path under different weights of the running cost. Moreover, the vehicle dynamics are manifested to have symmetry when the terminal pace is restricted to be equal to the initial one. We showed that the proposed model gives reasonable trajectories, and calibrated and validated it against real-world trajectory data.
Please note that simplified vehicle constraints on its motion are considered in our work. One may argue that the simple constraints for acceleration are not realistic. It can be improved in future work by incorporating powertrain limits and road adhesion coefficient, which are mainly relevant for automated driving. Moreover, the calibration and validation for the proposed vehicular movement model (for one vehicle) are conducted in this study. However, the interaction effects are left as hypothetical thought-experiments. In future research, the distributions of the parameters under various conditions should be calibrated in detail by using numerous empirical trajectory data. We anticipate that a realistic model to describe the vehicle dynamics at intersections can be developed on the basis of the proposed model, which can be used to assess and optimize intersection design schemes and signal control strategies.