text:
Benders' decomposition for freeway network design under endogenous autonomous vehicles demand

Autonomous vehicles (AVs) have the potential to provide cost-effective mobility options along with overall system-level benefits in terms of congestion and vehicular emissions. Additional resource allocation at the network level, such as AV-exclusive lanes, can further foster the usage of AVs rendering this mode of travel more attractive than legacy vehicles (LV). In this study, we propose an integrated mixed-integer programming framework for optimal AV-exclusive lane design on a freeway network which accounts for commuters' demand split among AVs and LVs via a logit model. We incorporate the link transmission model (LTM) as the underlying traffic flow model. The LTM is modified to integrate two vehicle classes namely, LVs and AVs with a lane-based approach. The presence of binary variables to represent lane design and the logit model for endogenous demand estimation results in a nonconvex mixed-integer nonlinear program (MINLP) formulation. We propose a Benders' decomposition approach to tackle this challenging optimization problem. Our approach iteratively explores possible lane designs in the Benders' master problem and, at each iteration, solves a sequence of system-optimum dynamic traffic assignment (SODTA) problems which is shown to converge to fixed-points representative of logit-compatible demand splits. Further, we prove that the proposed solution method converges to a local optima of the nonconvex problem and identify under which conditions this local optima is a global solution. The proposed approach is implemented on two hypothetical freeway networks with single and multiple origins and destinations. Our numerical results reveal that the optimal lane design of freeway network is non-trivial and can inform on the value of accounting for endogenous demand in the proposed freeway network design problem.

Introduction
As mobility demand increases, transport infrastructures remain under-utilized majorly due to various human factors e.g., slow reflexes, poor and heterogeneous driving behaviour, safety concerns etc. Automation in transport sector is gaining more attention than ever to minimize these human inputs to exploit transport resources in an efficient and sustainable way. This growing attention in automation translates to a projected growth of US$173B in global autonomous driving market by 2030 [30]. Potential benefits of autonomous vehicles (AV) include improved throughput [29,35,12], traffic safety [5], travel speed, energy consumption [20] and vehicular emissions [7].
Regardless of the benefits of AVs, one would be too naive to assume that AVs will be immediately adopted by legacy vehicle (LV) owners in near future. A more reasonable assumption would be the existence of a transition period where interactions between LVs and AVs exist, leading to a gradual increment in market penetration of AVs. This transition period will be crucial as safety might be compromised in mixed operations of LVs and AVs, especially in case of arterial networks, involving pedestrians, cyclists and signalized intersections. In these arterial networks with heterogeneous traffic, fewer AVs have been found to have a negative impact on the average travel speed and string stability of the traffic flow [34,31,22]. During this transition period, the safety concerns due to mixed vehicular interactions could be minimised with AV-exclusive lanes segregating AVs and LVs in a network, providing a smooth transition of the current transport system to an automation-heavy transport system. These dedicated lanes would be advantageous in the current situation of imperfect AV technologies, which pose critical problems in traffic safety during lane changing in a congested network with mixed traffic. At network intersections, reservation-based models have been shown to have the potential to increase intersection capacity [3,4,25,14]. More recently, Rey and Levin [27] proposed a hybrid network control policy in such networks with dedicated lanes that provide access to "blue phases" during which only AVs can traverse traffic intersections.
Yu et al. [39] adopted a microscopic traffic simulation method to investigate the efficiency and safety of mixed traffic on highways with AV-exclusive lanes. Although, the safety of mixed traffic was found to be worsened with low market penetration of AVs, the simulation results from the car-following models of this study showed an increment of up to 84% in throughput of the traffic network due to presence of AV-exclusive lanes. These dedicated lanes may also facilitate cooperative adaptive cruise control (CACC) contributing to better traffic flow performance [34], improved highway capacity [21], decrease in fuel consumption and emissions [24]. On freeways, CACC was found to significantly increase capacity with a moderate to high market penetration rate [29]. However, with a nave deployment strategy of AV-exclusive lanes, CACC was found to increase the total system travel time on freeways [19]. Talebpour et al. [32] simulated traffic flow under different penetrations rates of AVs on a two-lane and a four-lane freeway segment in Chicago, Illinois with mandatory and optional usage of AV-exclusive lanes. At market penetration rates of more than 50% for the two-lane highway and 30% for the four-lane highway, a potential benefit in terms of throughput and travel time reliability was observed in this study, with optional use of the AV-exclusive lanes yielding the most benefit. In another study, this optional use of AV-exclusive lanes was found to be beneficial only with a high market share of AVs [18].
Further, with the advent of AV technology in the market, it is critical to gain individual motivations for choosing to own or use AVs as a service. Hence, an endogenous demand model becomes imperative in such traffic flow modelling involving multiple vehicle classes. Based on a stated preference survey from 721 individuals living across Israel and North America, 44% of the sample population were found to prefer LVs over AVs [9]. Another stated choice survey for the adoption of shared AVs suggested service attributes including travel cost, travel time and waiting time to be critical determinants of the use of shared AVs [10].
To allow the transportation system to adapt with such disruptive technologies, we explore the potential benefits of introducing AV-exclusive lanes on freeways where the conflict points are considerably lesser than in local arterial networks and require minimal technological interventions. From previous studies, it is apparent that deploying AV-exclusive lanes in all links in a network may not be advantageous as it may significantly increase congestion and reduce throughput in the regular lanes; yet no systematic approach to identify the optimal location of AV-exclusive lanes has been proposed. In this work, we propose to investigate the problem of locating AV-exclusive lanes in a freeway network so as to reap maximum benefits of the deployed infrastructure.
We build on the literature of system optimum dynamic traffic assignment (SODTA) formulations to represent traffic dynamics, and obtain an analytical model amenable for exact optimization. Although a few studies [17,16,11,1] proposed SODTA approaches to model user route choice in general networks, this behavioral assumption is limiting, especially in the context of endogenous travel demand. Alternative approaches such as cross-nested and user-equilibrium models have also been proposed for characterizing different route choice principles [36,19], these research efforts are focused on static traffic assignment formulations or do not consider endogenous travel demand. To circumvent such assumptions, we focus on the design of freeway networks with single path per OD pair and we leave the study of general networks for future research. In this context, there are no route choice decisions, and the SODTA component of the formulation provides an analytical representation of traffic dynamics in the freeway network which maintains first-in-first-out (FIFO) conditions. Further, we do not model lane-changing behaviour of vehicles and do not account for vehicle holding in the network which is defined as the reluctance of vehicles moving forward from upstream to downstream links even with an availability of vacant spaces downstream [16,17].
For modelling the traffic flow in the network, we adopt the link transmission model (LTM) [38] and modify it to account for two vehicle classes: LVs and AVs. The LTM has been found to be useful to scale-up SODTA formulations for network design [1] as well as for solving routing problem for shared AVs in congested networks [11]. In the proposed formulation, we avoid mixed vehicular interactions of LVs and AVs by providing AV-exclusive lanes allowing AVs to behave autonomously only on AV-exclusive lanes and assume both LVs and AVs to behave identically on regular lanes. We capture the demand split of each vehicle class with a logit model, embedded in the proposed formulation. We model lane design decisions using binary variables, thus the proposed formulation results in a challenging non-convex MINLP.
We adopt a Benders' decomposition approach to implement this MINLP on a multi-OD freeway network which embeds a fixed-point algorithm in its sub-problem to account for endogenous travel demand.
This study makes the following contributions. We propose an optimization framework, involving an optimal design of dedicated lanes for AVs on a freeway network based on the demand of LVs and AVs which are endogenously calculated with a demand model embedded in the formulation. To the best of our knowledge, an integrated framework of this kind has not been explored in the literature. The computational complexity of such a nonconvex formulation is dealt with a Benders' decomposition approach in which the sub-problem embeds a fixed-point procedure to account for the endogenous demand model, retaining the linear structure of the underlying traffic flow model. We prove the existence of a fixed-point and our proposed algorithm always converges to this fixed-point. We also identify the conditions under which our algorithm is globally optimal. We implement the proposed algorithm on a freeway network which reveals that deploying a maximum of exclusive autonomous vehicles does not necessarily minimize network travel time.
This paper is organized into five sections. Section 2 presents the problem formulation followed by the solution methodology in Section 3. The proposed formulation is studied on two numerical networks, presented in Section 4. Section 5 presents the key findings of the study along with future research directions.

Freeway network design problem
In this study, we develop an MINLP model for optimal AV-exclusive lane design along with endogenous estimation of AV demand. We develop this MINLP in three stages. To begin with, we propose a linear programming (LP) framework for a lane-based LTM formulation with fixed AV-exclusive lanes, fixed proportion of AVs in the network and a system level objective. Further, we introduce a binary variable to obtain optimal lane design for AVs for an improved system performance, resulting in a mixed-integer linear program (MILP). Finally, we bring in a logit model to estimate the endogenous demand for each vehicle class which introduces non-linearity in the model, resulting in an MINLP. We circumvent this non-linearity by proposing a fixed-point algorithm along with method of successive averages (MSA) to obtain convergence of the fixed-point as explained in Section 3.
Let G = (N, A) be a directed network where N is the set of nodes and A is the set of arcs. In the proposed formulation, each node of the network is represented by a set of incoming and outgoing arcs. The set of arcs is partitioned into three subsets: source centroid connectors denoted as A r , sink centroid connectors as A s , and physical links formed by the remaining links of the set; i.e., A\{A r ∪A s }. The set of origin-destination pairs is represented by K. Γ − (i) and Γ + (i) represent the set of predecessor and successor links of link i. Table  1 presents the rest of the notations of the proposed formulation.
We choose the total system travel time (TSTT) as the metric of system performance, consisting of travel times of two vehicle classes: LV and AV. The objective function in the TSTT lv and TSTT av can be obtained from the difference in cumulative inflows and outflows of each link, representing the number of LVs/AVs present on that link at each timestep and the number of timesteps they spend on that link. The underlying LTM provides these cumulative inflows and outflows as output as explained in Section 2.2. The total vehicular demand in the network is the sum of LVs and AVs in the network as follows.

Endogenous demand model
We denote the time-varying total vehicle demand at time t between each origin-destination pair in the network by D o,d (t) and assume it to be fixed in our formulation. D o,d (t) is presented as a sum of the demands of two vehicle classes, LVs and AVs, in Eq. (2).
Though, the total demand in the network is fixed, the demand corresponding to each vehicle class varies depending on the proportions of AVs (p o,d ) between each OD pair in the network and they are obtained from Eqs. (3a) and (3b).
In the proposed formulation, the demand for each vehicle class is endogenous to the proportions of AVs (p o,d ), which is obtained based on the attractiveness of the modes in the network. We adopt a logit model to quantify this attractiveness, as summarized in Eq. (4).
We assume that the utility of each mode depends only on the average travel times of all the vehicles of that mode between each OD pair. The average travel times of LV and AV are denoted by τ o,d lv and τ o,d av and the utility of each mode is obtained from Eqs. (5a) and (5b).
Here, β represents the total amount of utility gained while making a trip, whereas, β o,d τ lv and β o,d τav are the disutilities created per unit travel time by LVs and AVs respectively. We refer to a study by Wong et al. [37] to obtain the values of β o,d τ lv and β o,d τav . We modify the logit model presented in Eq. (4) with the difference in utilities between the modes where the utility functions consist of the coefficients (β o,d τ lv , β o,d τav ) of average travel times of corresponding modes, as presented in Eq. (6).
The OD average link travel times of each vehicle class are obtained from the lane-based formulation of LTM which is the underlying traffic flow model in the proposed formulation.

Network dynamics
The LTM, proposed by Yperman et al. [38], is a numerical solution method for dynamic network loading, developed based on the first order kinematic wave theory [15,28]. In this study, LTM is chosen as it involves fewer variables per link compared to models such as cell transmission model (CTM) [2]. We adapt the conventional LTM to accommodate two vehicle classes by introducing two types of lanes: AV-exclusive and regular lanes. An AV-exclusive lane differs from a regular lane in terms of following headway, capacity and speed of backward shockwave propagation. The difference in these traffic flow characteristics affects the fundamental diagrams of traffic flow significantly as explained in the following subsection.

Fundamental diagram
The fundamental diagram of traffic flow reflects the relationship among the macroscopic traffic flow parameters of a network: traffic flow, density and speed. These relationships approximate all possible stationary traffic states during the analysis period and provide significant insights regarding the overall behaviour of traffic in a network.
The shape of the fundamental diagrams depicting these relationships may vary depending on the assumptions and approximations of a study. Greenshields et al. [8] was the first to propose a parabolic relationship between traffic flow and density. Later on, Newell [23] provided a simplified approach to the kinematic wave theory of traffic flow and developed a triangular shaped fundamental diagram, defined by three parameters: a fixed free-flow speed (v f ), capacity or maximum flow (q) and a jam density (K jam ). Yperman et al. [38] adopted this simplified fundamental diagram while developing the LTM which is the underlying traffic flow model in our formulation.
In a network with AV-exclusive lanes, the macroscopic traffic flow parameters may be significantly affected by faster reaction times of AVs leading to reduced following headway, increased throughput and faster propagation of backward shockwave due to congestion. Levin and Boyles [13] found considerable difference in the shape of the fundamental diagram for different reaction times of a characteristic vehicle. This study also showed how capacity and wave speed increase as the AV proportion increases with the human drivers having double the reaction time of AVs. Tientrakool et al. [33] demonstrated that due to tighter time and space headways among vehicles, the capacity of a lane could be approximately tripled by converting it into an AV-exclusive lane. Hence, while comparing traffic flow on AV-exclusive and regular lane, the shape of the triangular fundamental diagram will be significantly different due to the changes in capacity (q) and backward wave speed (w) leading to the same jam density (K jam ).
In this study, we consider the AV-exclusive lanes to have double the capacity (q av = 2q lv ) and backward wave speed (w av = 2w lv ) of that of the regular lane while keeping the free-flow speed (v f ) and jam density (K jam ) the same for both lane types. The fundamental diagrams of traffic flow on both of these lane types is shown in Figure 1.

Traffic flow propagation
The LTM keeps track of the vehicular flow in the network with cumulative inflows and outflows of each link at each time-step. We develop a lane-based LTM where the vehicle class-specific cumulative inflows and outflows from each lane i towards destination k at time t are denoted by z k+ i,lv (t) z k+ i,av (t) and z k− i,lv (t) z k− i,av (t) respectively. The demand corresponding to each vehicle class is loaded into the network as the cumulative inflow to the source centroid connectors as shown in Eqs. (7a) and (7b).
The cumulative inflow to the other lanes at time t is defined as the sum of transfer flows from all the incoming lanes predecessor (Γ − ) to that lane over all the timesteps up until t. We define these in Eqs. (8a) and (8b).
Similarly, the cumulative outflows (z k− i,lv (t), z k− i,av (t)) from a lane at time t is defined as the sum of transfer flows to all the outgoing lanes successor (Γ + ) to that lane over all the timesteps up until t. We define these in Eqs. (9a) and (9b).
The LTM has been built based on three flow components: sending flow, receiving flow and transfer flow. Sending flow is defined as the amount of vehicular flow allowed to go out from link i to link j respecting its flow capacity. Yperman et al. [38] Eqs. (11a) and (11b) present the capacity constraints on sending flow for regular and candidate AV lanes respectively. On regular lanes, the total flow of LVs and AVs is restricted to the capacity of regular lanes (δq lv ) whereas, on candidate AV lanes, we introduce a binary In the LTM, the receiving flow is defined as the amount of vehicular flow allowed to be received at link j from link i depending on the congestion level and the capacity of link j.
The receiving flow constraint, as presented in Eq. (12), is derived based on the backward propagation of a congested traffic state from the downstream boundary of a link which reaches the upstream boundary L i wav time units later. Here, w av denotes the backward wave speed of the congested traffic state.
Similar to Eqs. (11a) and (11b), Eqs. (13a) and (13b) represent the capacity constraint on receiving flow of a link with the binary parameter, b, depending on the lane being a regular or candidate AV lane.
In the proposed formulation, LVs are restricted from entering an AV-exclusive lane. We formulate integer-linear constraints to implement this restriction in our model as presented in Eqs. (14a) and (14b). Using the binary lane design variable (b i ), the transfer flow of LVs at any timestep is either restricted or kept free, i.e., equal to the capacity of regular lane, for a downstream AV-exclusive or regular lane respectively.
Eqs. (15a) and (15b) conclude the lane-based LTM formulation ensuring the exit of all the vehicles that entered into the network and reaching their respective destinations at the end of the last timestep (t).

MINLP formulation
In the LTM, the cumulative inflows and outflows of each lane at each timestep track the vehicular flow in the network. The difference between these inflows and outflows of a lane at each timestep represents the number of vehicles present in that lane for δ time units, where δ is the duration of each timestep. Hence, the sum of these differences over all the lanes, OD pairs and timesteps will provide the TSTT of each vehicle class in the network as follows.
The average travel times are estimated based on the number of timesteps each vehicle spends on each link, averaged over the total demand of that vehicle class as presented in Eqs.
We rewrite the objective function of the proposed MINLP presented in Eq. (1) as follows.
Note that, as LVs are restricted on AV-exclusive lanes, we fix at least one path with regular lanes between each OD-pair in our model for movement of LVs. The resulting MINLP formulation FNDP represents the freeway network design problem. The remaining variables are used to account for congestion effects and endogenous travel demand. Travel demand is loaded into the network through the source centroid connectors as expressed in Eqs. (3a) and (3b). At the completion of the trips, Eqs. (15a) and (15b) ensures that vehicles leave the network through the sink centroid connectors. On a freeway, these source and sink centroid connectors represent on-and off-ramps respectively. If the network is unable to accept demand due to congestion on the freeway, vehicles may be held at on-ramps which are assumed to have sufficiently large capacities. Since waiting time is penalized in the objective function, the outputs of the proposed formulation can be interpreted as the level of control at the freeway on-ramps.
We analyse the proposed model in the following sections.

Fixed-point analysis
To motivate the design of a dedicated solution method and to provide insights into the behavior of FNDP, we consider a simplified version of the model wherein the endogenous demand p and the lane design b are assumed fixed. This simplified model is called SP(b, p) and presented below.
The variables involved in Model SP(b, p) are transfer flows y ∈ Y and cumulative inflows and outflows z ∈ Z. We consider a single-OD network with two links with fixed AV-exclusive lanes as illustrated in Figure 2, where the fixed AV lanes are shown in blue. The first link consists of three regular lanes and one AV-exclusive lane followed by a capacity drop on the second link which has one LV and one AV-exclusive lane. We solve SP(b, p) for a series of OD proportions of AVs (p o,d ) and calculate the corresponding logit-derived proportion of AVs (p o,d logit ) using Eq. (6) based on the optimal solution of SP(b, p). Note that if p = p logit , then the demand splits are logit-compatible, i.e. equilibrated, and this solution corresponds to a fixed-point. We define this fixed-point as follows.
Definition 1 (Fixed-point). Let F : P → P be a continuous function of the OD proportion vector p ∈ P. We say that p is a fixed-point if F (p) = p. The regular and AV-exclusive lanes differ from each other in terms of capacity and backward wave speed of congestion propagation. In this study network, the capacity (4320 veh/hr) and backward wave speed (28.4 km/hr) of an AV-exclusive lane is taken as double that of an regular lane due to the inter-connectivity of AVs leading to better traffic flow and faster congestion propagation. The length and jam density of the links are 800m and 200 veh/km respectively with a free-flow speed of 90 km/hr for both vehicle classes. The demand is loaded through the source centroid connector 1 into the network. The capacity and jam density of the source and sink centroid connectors are set to very high values with a negligible length for instantaneous loading of demand into the network based on available network capacity. The values of these network parameters are provided in Table 2.    We observe that all fixed-points have equal TSTT. This experiment highlights that fixed-points may be non-unique, but may also correspond to identical TSTT. From Figure 3, we observe that in each case, the line plot depicting the relationship between p and p logit crosses the reference line at least once, referring to the existence of at least one fixed-point satisfying the logit model. We prove this existence of at least one fixed-point with Proposition 1. Proof. We show that we can construct such a continuous function F (p). Let h(p) be the function mapping the OD proportion vector p to the optimal solution of the linear program SP(b, p) for fixed lane design vector b as defined by (2). Let z ∈ Z be the vector of cumulative inflow and outflow and let y ∈ Y be the vector of transfer flows obtained after solving SP(b, p), formally: Let g(z, y) be the function mapping the optimal cumulative inflows and outflows, and transfer flows to average, class-based (LV and AV) OD travel times (τ lv , τ av ) ∈ T as defined in Eqs. (17a) and (17b): Finally, let f (τ lv , τ av ) be the function mapping average class-based OD travel times to OD proportion p ∈ P via the proposed logit model as defined in (6): Let F (p) = f (g(h(p))), F is a continuous function from the compact convex set P to itself. By Brouwer's theorem, there exists at least one fixed-point p such that F (p) = p.
In Section 3, we develop the solution methodology to solve the non-linear Model 1 with variable lane design and endogenous demand based on the embedded logit model.

Benders' decomposition approach
In this section, we propose a decomposition approach with variable lane design and endogenous demand for each vehicle class. The purpose of this development is to identify the crucial locations in a network where providing an AV-exclusive lane will reap the maximum benefit for a given AV demand. As LVs are restricted on AV-exclusive lanes, it is necessary to deploy AV-exclusive lanes judiciously to cater to both vehicle classes keeping the social welfare into perspective. We introduce Benders' decomposition method for this purpose which iteratively explores possible lane designs in a master problem and, at each iteration, solves a sequence of SODTA problems which is shown to converge to fixed-points representative of logit-compatible demand splits.
Benders' decomposition approach eases up the computation burden of a mathematical model by partitioning the overall formulation into a relaxed master problem with mainly integer variables and a subproblem with all the continuous variables. For a detailed review of Benders' approach, one can refer to Rahmaniani et al. [26]. For problems with minimizing objective function, such as the proposed model in this study, the relaxed master problem provides a lower bound at each iteration of Benders' method. Whereas, the subproblem handles the complicated constraints of the original problem which is solved iteratively for each relaxed solution of the master problem. At each iteration until convergence, the Benders' method generates either a feasibility cut or an optimality cut towards obtaining the optimal solution. A feasibility cut is generated to eliminate an infeasible solution provided by the subproblem, preventing the master to produce it again. On the other hand, if the sub-problem solution is feasible, an upper bound is obtained and an optimality cut is derived towards closing the optimality gap.
In our decomposition of Model FNDP, we only retain binary lane design variables b in the master problem MP, which is summarized below: In the proposed model, the subproblem (SP(b, p)) is initiated by fixing the set of binary lane design variables. Depending on the feasibility of the subproblem, dual prices or rays for each constraint are calculated, followed by solving the master problem which provides the values for the next set of binary lane design variables. This iterative process continues until we reach an exact solution of the objective function.
The subproblem in the proposed formulation consists of two components: the lanebased LTM with system-level objective function and an endogenous demand model. As the endogenous demand model introduces non-linearity in the formulation, we develop a fixed-point algorithm for solving the subproblem as explained in the following subsection.

Fixed-point algorithm
The endogenous demand model is crucial to study the effect of infrastructural changes such as AV-exclusive lanes on AV demand. Note that, the formulation without the endogenous demand model is a useful model by itself as it can estimate the progressive deployment of AV-exclusive lanes in a network corresponding to incremental penetration of AV demand. However, this model does not account for the effect of these AV-exclusive lanes on the demand of each vehicle class.
We adopt the logit model for endogenously estimate the demand for each vehicle class. As the logit model is nonlinear, we use a fixed-point algorithm to separate the endogenous demand component from the SODTA formulation. The fixed-point algorithm estimates the proportion of AVs (p o,d ) in an iterative process involving the logit model, presented in Eq. (6). The Algorithm 1 presented in the next subsection, keeps the nonlinear logit model out of the MILP, keeping it a linear mathematical formulation. The new proportion of AV demand is obtained based on the difference in utilities between the modes, substituted in the logit model (as shown in Eq. (6)), which is fed back to Model 2 for subsequent iterations. We adopt the method of successive averages (MSA) for convergence of the fixed-point algorithm which is based on a predetermined move size along the descent direction. Note that MSA is guaranteed to converge after a certain number of iterations by providing lesser weightage to subsequent solutions at each iteration. This iterative process of MSA may disregard the instability in the fixed-point solution. We monitor this instability by checking the value of p o,d before and after implementing MSA.
The Benders' decomposition algorithm along with fixed-point algorithm is presented in Algorithm 1. We next prove that Algorithm 1 converges to a local optima of Model FNDP in Proposition 2 and identify under which conditions this local optima is a global solution in Proposition 3.  Proof. If all fixed-points have equal TSTT, then Benders' cuts are guaranteed to never overestimate the lower bound Z obtained from solving the master problem MP. Hence, all local optima of Model FNDP have equal and minimal TSTT.

Algorithm 1: Benders decomposition with fixed-point algorithm for
(Solve logit model with free-flow travel times) Proposition 3 identifies the conditions under which the solution of the proposed algorithm is globally optimal. Although we are not able to mathematically prove that in the presence of multiple fixed points all fixed points yield identical TSTT, in our investigations we consistently observe this behavior on congested freeway networks.
In the next section, we implement this algorithm on a single-OD and a multi-OD network.

Numerical experiments
The proposed formulation is implemented on the same single OD network ( Figure 2) presented in Section 2 along with a multi-OD freeway network ( Figure 6).

Single OD freeway network with fixed lane design
We implement the solution methodology presented in Section 3 on the single OD freeway network with fixed lane design for AVs. As mentioned earlier, the presence of the binary lane design variable makes the proposed formulation a mixed-integer linear program, whereas, the endogenous demand model introduces non-linearity due to the structure of the logit model. We develop the fixed-point algorithm to circumvent non-linearity and Benders decomposition method to handle the binary lane design variable in the master problem. To understand the performance of the fixed-point algorithm alone, we first implement Model 2 on the single OD network with fixed dedicated lanes for AVs. We initialize the algorithm with a proportion of AV obtained from Eq. (6) substituting the free-flow travel time between the OD pairs. Model 2 is solved at each iteration of this algorithm. The value of p o,d is updated at each iteration with the proportion (p logit ) obtained by the endogenous demand model until convergence where a fixed-point is reached (p = p logit ). The network characteristics of this test network is same as presented in Table 2. A time-varying demand profile is selected for this analysis which is loaded into the network every 2 minutes for the first 8 minutes of a total analysis period of 100 minutes. The total demand is 2950 vehicles which includes both LVs and AVs. Note that, the objective function of the proposed formulation takes into account the waiting time of vehicles at the source centroid connector depending on the available capacity in the network. We consider a timestep (δ) of 30 seconds at which the cumulative inflows and outflows of each link are updated by the underlying LTM. For a single OD network with fixed dedicated lanes for AVs with β o,d τav = 0.0018, we plot this convergence of the fixed-point algorithm in Figure 5. Figure 5 shows that the algorithm converges to a fixed-point after 67 iterations at p = 0.68 which is one of the fixed points presented by the purple dot in Figure 4. At convergence, the value of the objective function (TSTT) is found to be 2604430 veh-sec with a computation time of 24.067 seconds.
In this case study on a single OD network with fixed AV-exclusive lanes, we observe that the fixed-point algorithm performs well with fast convergence. In Subsection 4.2, we introduce the Benders' decomposition method to handle variable lane design problem on a multi-OD network.

Multi-OD freeway network
We implement the proposed model with Benders method and fixed-point algorithm on a multi-OD freeway network, presented in Figure 6. This 27km long freeway consists of 6 OD pairs, 9 links, 4 on-ramps and 4 off-ramps. Each of these links has either 2 or 3 lanes as shown in Figure 6 where one lane in each link, showed in green colour, belongs to the set of candidate AV lanes (A av ) which could potentially converted into an AV-exclusive lane for better system performance. The total vehicular demand is considered to be 3000  We design 7 experiments to study the performance of our algorithm. Apart from the base case, the proposed algorithm is implemented for different demand (±25%), capacity of AV lanes (±25%) and coefficient of average AV travel times (±25%). The performance of Algorithm 1 for these different scenarios is presented in Table 3. We observe that it is not beneficial to convert all candidate AV lanes to AV-exclusive lanes. For example, 7 out of 9 candidate lanes are converted to AV-exclusive lanes for the base case which remains the same for the reduced demand case as well. However, for an increased demand, a lesser number of lanes are converted to AV lanes (6). This could be due to increased LV demand which requires more regular lanes in the network. In case of an increased AV lane capacity, the system performance is significantly enhanced by allocating more AV-exclusive lanes (8) compared to the reduced AV-lane capacity case (5). On the other hand, the coefficient of average AV travel times is found to have no effect on AV lane design in the network. These deployed AV-exclusive lanes are found to decrease the TSTT by around 10% in all the scenarios, with a maximum improvement of 12.71% for increased AV-lane capacity. Hence, AV-exclusive lanes are always found to have a positive impact in network performance.
The computation time of the proposed algorithm varied from around 0.5 to 4 hours on an Intel(R) Core(TM) i7-6700 @3.40GHz CPU with 16GB RAM. The base case was found to converge within 2 hours with 55 iterations of Benders and a total of 269 fixed-point iterations. Interestingly, the change in the capacity of AV-exclusive lanes is found to have the maximum effect on computation time of the proposed algorithm with a 10 times increment in computation time for inflated capacity compared to the deflated capacity of AV lanes. Similar increment in the number of Benders and fixed-point iterations is also observed in these cases with inflated and deflated capacities. Table 4 summarizes the optimal proportion of AVs (p o,d ) obtained from the endogenous demand model in different scenarios. We initialize the values of p o,d by substituting freeflow travel times of vehicles in the logit model Eq. (4). Due to the structure of the logit model, the initial p o,d values are found to be inversely proportional to the free-flow travel time between OD pairs and vary from 0.28 to 0.44. As the free-flow travel times on regular and AV-exclusive lanes are the same, this variation is solely due to the difference in the coefficient of travel times in the logit model. In congested traffic conditions the proposed algorithm estimates these p o,d to be around 0.57 to 0.83. In the base case, the OD pair (3,17) is found to have the maximum value of p o,d with a value of 0.76, whereas the OD pair (3,9) requires the least proportion with a value of 0.63. With a reduction in total demand, this proportion is found to be reduced by a maximum of 10%. Whereas, with increased demand, the maximum increment in p o,d is only around 4.6%. This could be due to the increased total demand which require fewer AV-exclusive lanes, as shown in Table 3, to cater for the increased LV demand. Hence, increasing p o,d may not provide much improvement in system level performance. On the other hand, p o,d values are not found to be affected significantly with inflated and deflated capacity of AV-exclusive lanes. However, the different coefficients of the average travel times of AVs (β av ) are found to considerably affect the AV proportions. A reduction in this coefficient means that average travel time will create less disutility for AVs rendering this mode more attractive than AVs and vice-versa. The output from the proposed algorithm shows similar effect with a maximum increment of 12% in p o,d for a 25% reduction in β av and a maximum of 15.7% decrease with an increment in β av . We present the average travel times LVs and AVs on each link for all the 7 experiments in Tables 5 and 6. As mentioned earlier, the proposed algorithm triggers 7 out 9 candidate AV lanes into AV-exclusive lanes. These exclusive lanes are found to have a significant effect on the average travel times of LVs with a maximum reduction of 45% for the base case. Similar improvement is observed for other demand scenarios as well with a maximum reduction of 41 and 46% respectively for LVs. These maximum reductions are mainly observed on the waiting times of vehicles on on-ramps. With a deflated capacity of the proposed AV lanes, the average reduction of average travel times on congested links is found to be around 8% with a maximum reduction of 23% in the first on-ramp of the network. Whereas, the inflated capacity provided a average reduction of 12% in the travel times with a maximum reduction of 52%. In the case with reduced β o,d τav value, the average reduction of travel times for LVs are found to be 5.6% compared to the case with increased β o,d τav value where this reduction is 14%. This output is quite intuitive as the reduced a β o,d τav will create less impedance on the mode choice decisions of users based on AV travel times, leading to more AVs in the network and less travel time benefits of for LV users.
With the proposed AV-exclusive lanes, the average travel times of AVs are found to be reduced by an average of 15%. Note that, these average travel times correspond to AVs on both regular and exclusive lanes. The maximum reduction of 60% is observed for the waiting time on source connector 10. The average reductions in travel times are around 16.4% and 13.4% for the reduced and increased demand cases respectively. The deflated and inflated capacities of the AV-exclusive lanes reduced the average AV travel times of AVs by 12.8 and 120 120 120 120 120 120 120 120 120 120 12 120 120 120 120 120 120 120 120 120 120 14 120 120 120 120 120 120 120 120 120 120 16 120 120 120 120 120 120 120 120 120 120 17.5% respectively with a maximum reduction of 64% for the inflated capacity case. With an decreased and increased β o,d τav the average reduction in the average travel times are around 16.4 and 15% respectively.  120  120  120  120  120  120  120  120  120  14  120  120  120  120  120  120  120  120  120  120  16 120 120 120 120 120 120 120 120 120 120

Conclusion
In this study, we proposed an optimization framework to solve a multi-OD freeway network design problem for optimal lane design under endogenous AV demand. Due to the presence of binary lane design variables and the endogenous demand model, the proposed formulation results in a nonconvex MINLP. We tackle this challenging problem by introducing Benders' decomposition approach which iteratively explores possible lane designs in the master problem and at each iteration solves a sequence of SODTA problems which is shown to converge to fixed-points representative of logit-compatible demand splits. The proposed approach is implemented on two hypothetical freeway networks with single and multiple origins and destinations. On the single-OD network, the fixed-point algorithm is found to converge at multiple fixed points providing different proportions of AV demand. However, these multiple fixed-points are found to have no effect on the objective function of the problem. Here, we prove that for a fixed lane design, there exists at least one fixed-point representing the proportion of AV demand in the network. We also prove that the proposed solution method converges to a local optima of the nonconvex problem and identify under which conditions this local optima is a global solution.The numerical results on the multi-OD network show that it is not beneficial in terms of system performance to provide AV lanes for all the links in the network.
We observe that the optimal lane design of freeway network is non-trivial and can inform on the value of accounting for endogenous demand in the proposed freeway network design. The proposed model may also be useful for designing ramp metering for a freeway network with AV-exclusive lanes. The following limitations may be examined in future research. The proposed algorithm is designed for a multi-OD freeway network with single path between each OD pair, hence route choice modelling is averted. We do not model lane-changing behaviour nor attempt to model vehicle holding issues. This work can also be extended by investigating user equilibrium (UE) vs SO route choice in an arterial network with dedicated AV lanes incorporating mixed vehicular interactions.


An image vision and automatic calibration system for universal robots

The rapid development of technology is causing the replacement of many traditional manufacture industries by automation. Robotic arms are now commonly used in many sectors. The requirements for robotic arms in different sectors are quite different; but, in general, robotic arms not only save cost and man power, they also improve safety. The aim of this study was an investigation of the integration of image identification with robotic arms. To do this, the Denavit–Hartenberg transformation matrix was used to analyze the mechanical kinematics of the joints of each robotic arm axis. This allowed the spatial relationships between the Cartesian 3 D coordinate system and the joint of each axis to be determined and communication between the robotic arm and image identification to be established. A robotic arm prototype platform with automatic image identification and calibration was constructed using a quick and robust method. Several of the variables that exist in real robotic arm applications have been solved in this study: primarily the accuracy error that occurs when repeated gripping of workpieces is done and a movement and placement track is set up. Deviations occur frequently and if image identification can be applied to offset repeated inaccuracy, the quality of finished products will be more consistent. The Universal Robots 5 six-axis robot and cameras, which use a filter, Sobel computation, and Hough transformation computer vision image processing technology, together comprise a system that can automatically identify the workpiece and carry out loading and unloading accurately.

Introduction
The purpose of the use of robotic arm is to effectively reproduce the functions of a human arm and to carry out operations under automatic control. These arms are now commonly used on many different industrial production lines. Their primary function is to carry out repetitive work more efficiently and faster than humans. 1,2 The most obvious difference between a robotic arm and normal automatic machinery is the multiple links and joints. With these functions, the robotic arm can freely move in three-dimensional (3D) space. Most robotic arms have four main parts: the main machine, controller, servo mechanisms, and sensors. The robotic arm used in this paper was a UR5 six-axis robot (Universal Robot 5). 3,4 The six-axis control involves the kinamatics 5 of the Denavit-Hartenberg (D-H) transformation matrix, 6,7 dynamics, as well as control theory, 8,9 including the computation of six-axis rotational angles from the work point and pose of the robotic arm and the derivation of spatial work point position and pose from rotational angles to complete spatial track planning, 10 etc.
The Hough transformation was used for image feature detection and identification and is widely used in computer vision, graphic analysis, and digital image processing. [11][12][13] Its algorism process is roughly as follows: when an object needs to be identified, the type can be found from the image and easily derived using a standard Hough transformation.

Robotic arm kinematics
The first step in robotic arm operation is accurate marking of the target object position as well as the relative position of the robotic arm to the target object. In other words, the orientation relationship between them must be determined. After gaining the necessary information about the coordinates, each stage of coordinate transformation can be carried out. Kinematics is used to describe the tracking of an object in motion and involves forward and inverse kinematics. 17 Forward kinematics determines the robotic arm pose and end point via the rotational angle of each axis joint. Inverse kinematics back calculates the joint angle of each axis from the robotic arm pose and end point position. Forward kinematics is the basis of inverse kinematics and can be used for checking the correctness of inverse kinematics. In industrial applications, a series of end point positions are often given. Inverse kinematics is used to determine the corresponding joint angle of each axis to each position and to control the motors, so that each joint can reach the determined position. This allows the robotic arm to complete each given pose. The kinematic relationship is shown in Figure 2. The D-H forward kinematics robotic arm coordinate system Forward kinematics involves sending position angle information to the controller for each joint to establish a position. In the D-H coordinate method, appropriate parameters are defined which make originally complex transformation simple and clear. The UR5 robot positions and local reference coordinates of each axis joint [18][19][20][21] are shown in Figure 3.
Standard D-H parameter, i-robotic arm axis joint, h-joint angle, d-linkage offset, a-linkage length, a-linkage torsion, unit: mm, or rad are shown in Table 1. 8,18,21,22 Via the D-H transformation matrix, equations shown in equations (1) to (4) Therefore, the result of each axis via the coordinate system transformation matrix iÀ1 i T can be determined from the above four transformation matrices as shown in equation (5)   It is shown in equation (6) after rearranging equation (5) iÀ1 By combination with the (Table 1) UR robot D-H parameters, we can determine the transformation matrix of each axis joint as shown in equations (7) to (12) Finally, the transformation matrix is determined as shown in equation (13) After simplification, the parameters are as in equation (14) s where parameters are shown in equation ( When all h are known, from equation (13), we can obtain the UR5 robot rotational matrix R as in equation (16) and the final coordinate matrix P as in equation (17) R ¼ n x o x a x n y o y a y n z o z a z

Robotic arm inverse kinematics
If the pose of the robotic arm end point relative to the machine base coordinate is known, we can back calculate the joint angle h i of each axis, by inverse kinematics. 9 First, define P u v as the vector of the origin coordinate of the joint of the vth axis to that of the joint of the uth axis. Also, use the coordinate of the joint of the uth axis to represent its vector is shown in equation (18) Since P 0 6 and 0 6 T are known, we can determine P 0 5 as shown in equation (19) from Table 1 P Since P 0 6 , P 0 5 , and 0 6 T are known, Table 1 can be used again, as well as the auxiliary angles of u and w from Figure 4, to determine h 1 as shown in equations (20) to (22) 23,24 Since P 0 6 , P 0 5 , 0 6 T, 0 1 T, and h 1 are known, we can determine P 1 6 as shown in equation (23) from Table 1 P 1 6 ¼ 0 6 T À Á À1 P 0 6 (23) From Figure 5, we can determine h 5 as shown in equation (24)   Since P 0 6 , P 1 6 , P 0 5 , 0 6 T, 0 1 T, h 1 , and h 5 are known, Table 1 is used again to determine 1 6 T as shown in equation (25) 1 where a, b, and c are shown in equations (26) to (28) Therefore, when h 5 is neither zero nor p, we can derive equation (29) as When h 6 is known, 5 6 T can also be known. When 5 6 T and h 5 are known, we can determine 4 5 T, and equations (30) to (32) can be derived Now, P 1 3 , P 1 4 , P 0 6 , P 1 6 , P 0 5 , 1 4 T, 0 6 T, 4 6 T, 0 1 T, 4 5 T, 5 6 T, h 1 , h 5 , and h 6 are known. Using Figure 6, we can derive h 3 and h 2 as shown in equations (33) to (34), where e, d, and n are auxiliary angles When h 2 and h 3 are known, 1 2 T and 2 3 T will be known. Therefore, we can obtain equation (35) Last, using equations (36) and (37), we can determine h 4 as shown in equation (38)

Analysis of the kinematic singularity
The six-axis robotic arm is driven by six individual motors that provide axial motion at each axis. The arm has six degrees of freedom in 3D space. Theoretically, the end point position of the six-axis robotic arm can reach any position and angle in 3D space. However, it can still get stuck in certain positions because some singularities remain. When it is located at a singularity, inverse kinematics can fluctuate. This poses a potential danger to later tracking planning. When the end point is near a singularity, a minute displacement change can cause a dramatic angle change to some axes. This can generate an approximate infinite angular velocity. If the inverse kinematics of the arm has a solution in the analytic form, 24 then it is possible to find the possible singularities by judging the abnormal condition where the denominator of the algebraic form of the solution is zero. When the value of the matrix determinant of the arm is zero, the robotic arm is located at a singularity. If the second, third, and the fourth axis are in parallel, in an arm with six degrees of freedom, the solution in analytic form that exists in inverse kinematics is given by equation (29), and when s 5 ¼ sinh 5 ¼ 0, calculation will show the value of the robotic arm matrix determinant to be zero. This means that when sinh 5 ¼ 0, the robotic arm is positioned at a singularity. Now, h 6 can be any value to fulfill the given constraint. For this, it is necessary that processing of h 6 near the singularity be smooth to avoid fluctuation h 6 and also to guarantee motion accuracy. As mentioned, the singularity can be solved using fuzzy motion planning, 25 constrained quadratic optimization, 26 intermediate desired value approach, 27 and the evolution algorithm. 28

Hough transformation and mass center calculation
In this study, the Hough transformation method was used to find the location of the image of the workpiece to be gripped and also to determine the mass center of the workpiece. After processing, this information was transferred to the robotic arm controller to complete the gripping task. The Hough transformation, also called Hough gradient method, is a special detection method widely used in digital image processing and computer vision analysis. The Hough transformation uses binarization to identify specific lines in the image. The theory is the location of line parameters using the positions of all the points in the image. One-to-many mapping of each point in the image to the parameter space is used to get the possible parameter values. The parameter values generated by all points are accumulated to determine the line parameter that presents most obviously in the parameter space. The thinking behind the Hough circular transformation [29][30][31][32] is that each non-zero pixel point is likely to be a point on a circle, and an accumulated coordinate plane is generated by voting. The position of the circle is determined by setting up an accumulated weight. On the line coordinate, the circular equations of all points on the same circle are the same. Therefore, when mapping to the 3D space Cartesian coordinate system, also known as the abr-coordinate system, they will be the same point. Therefore, by judging the accumulated intersections of each point in the abr-coordinate system, the points greater than a certain threshold value are circles. Any straight line that passes and holds a right angle to the circumference will intersect and be perpendicular to a point. The intersection of the points is the center of the circle. The threshold value can be judged from the number of intersections on the circle center, as shown in Figure 7.
For a standard Hough transformation, 33 three parameters are needed to describe a circle, the circle center, and the radius (a, b, r). Any point (x, y) in the image space can be mapped to 3 D space (a, b, r) to form a circular cone as shown in equation (39) Look for a possible circle center for the point on each edge and build a 3D circular cone accumulator by using the edge as the circle center. These circles will overlap on the target circle center. Last, find the circle center and the radius from the point with most overlap. The circle can be expressed as in equation (40) ðx À aÞ 2 þ ðy À bÞ 2 ¼ r 2 Analyze any pixel point in the graph corresponding to the graph in 3D space. If the radius r is certain, the pixel point corresponding to the 3D space will be a circle. If the radius r is not certain, then any shape with the radius r must be a circle. This circle will form a 3D circular cone as shown in Figure 8.
If a line is drawn along the gradient direction from any point on the circle, then the intersection of all lines may possibly be the circle center. First, calculate the possible circle center. Then, calculate the possible radius with the circle center. The circle that is checked must have a boundary. Detect the edge with the Sobel operator. 34,35 Technically, the Sobel operator is a discrete difference operator used to find the approximate gradient value of the graphic brightness function. This operator includes two sets of 3 Â 3 matrices, which are lateral and longitudinal. The Hough circular transformation process has two steps. First, the results and procedures of circle detection are shown in Figures 9 and 10. Then, the procedures and results of radius detection are shown in Figures  11 and 12.
The original color image (Figure 9(a)) processes the original image into a grayscale image (Figure 9(b)). Then, the Sobel operator is used to get the horizontal and vertical gradients (see Figure 9(c) and vertical (Figure 9(d)). Finally, the horizontal and vertical gradient values are combined to find the center point (Figure 9(e)).      ð Þ accumulators are compared with the threshold to find the circle in the image. Figure 11 illustrates the Circle Radius and Circle center Estimation Flowchart. The radius of each gradient value obtained from the gradient image is detected, and the value obtained from detection is saved in the N r ð Þ radius accumulator and compared with the threshold value to obtain the most suitable radius and the center of the circle. Movement in the positive or negative directions of the gradient is made to carry out detection within the range of the radius, see Figure 12(a). The points recorded after detection are assumed to be the center of circles, see Figure 12(b), and comparisons will result in the strongest being left, see Figure 12(c). Last, the radius using some of the assumed circle centers is calculated to find the optimum circle (Figure 12(d)). Last, Hough transformation is used to obtain the space graph, see Figure 12(e).
After obtaining the optimum circle, we can get the planar transform relationship between the workpiece to be gripped and the current image. The center of this circle can also be defined as the mass center. The mass center is the reference point for control, and its 2D planar coordinates represent the position of the workpiece to be gripped. The robot arm can now be directed to complete the gripping task.

Communication method-socket
To satisfy the communication requirements 36,37 for the robotic arm system, and to achieve support for some TCP/ IP functions that are not easy to integrate into the prevailing function library (TCP/IP Protocol Suite), we define some system call functions-Socket Interface. The Socket web framework is shown in Figure 13. The TCP/IP Socket protocol used in this study forms a satisfactory link between the PC end and the controlled end of the robotic arm.
The application layer is similar to that used by other network applications such as Telnet, FTP, WWW, POP Email, News, and BBS. The transport layer is responsible for network connection, division, and combination for transmitting messages and for providing control for data traffic between usage nodes. It also decides the service quality and provides reliable and efficient connection for network application nodes. The IP network layer receives the packets sent from the transport layer. Transmission is automatically routed to the address, and packet format can be changed to that of a different protocol, the network traffic condition is monitored, and an overall network topology framework is dynamically built to provide the best route for data transmission. The network access layer receives the packets sent from the IP network layer to divide, combine, and detect the data frame in more detail and correct errors, and ACK is used to determine whether the data is transmitted normally and controls the transfer rate. It is responsible for transmitting the original network data formed between 0 and 1. The method used in this study was TCP/IP. Socket is a port used for network communication in the UNIX system. In general program operation, it can be regarded as an operation for File I/O. Basically, this relies on two key points: read and write. How to read the incoming data and how to write the output. The basic steps before performing these operations are shown in Figure 14.

Experimental results and discussion
Details of the experimental flowchart are shown in Figure 15. An image of the object is acquired by the CCD camera and transferred to the PC via the USB connection. The Hough transformation method is used to analyze and compute the information until the workpiece to be gripped is found as shown in Figure 16. The actual flow is from the original RGB image input (Figure 16(a)) which has been reduced to grayscale (Figure 16(b)) to reduce its computational complexity and the SOBEL level (Figure 16(c)), vertical detection of the image (Figure 16(d)) has been started and integrated (Figure 16(e)), the gradient values have been converted into the Hough space ( Figure 16(f)) to determine if it is a circle, and finally, the circle identification ( Figure 16(g)) has been completed. The mass point of the object is then found by calculation and the D-H coordinate system transfers this information (by the socket method) to the UR5 robot controller. The robotic arm carries out the necessary movement to grip the workpiece. It automatically moves to the correct position and completes the movement as shown in Figures 17 and 18.
After the workpiece has been gripped, the robot arm receives commands (via Socket) to move to the fixture position where the workpiece is to be placed as shown in Figure 19. The second CCD camera provides  information for visual identification, and any deviation between the workpiece being held by the arm and the fixture will automatically be corrected. After this, the work continues, and the workpiece is placed smoothly into the fixture. The workpiece is finally picked out of the fixture after machining has been done. The system described here also has a user-friendly man-machine interface for monitoring the operations. The Hough transformation is used to identify the workpiece to be clamped on the captured image ( Figure 17(a)). The robot arm moves to waiting point 1 (Figures 17(b) and 18(a)). After the centroid of the workpiece is found, the coordinate of the workpiece to be selected is sent to the UR5 robot arm using the Socket method, and UR5 robot arm is moved to the first position, point 1. After a wait at point 2 (Figures 17(c) and 18(b)), it finally moves to the workpiece position (Figures 17(d) and 18(c)) to capture the workpiece (Figures 17(e) and 18(d)).
The Hough transformation is useful for recognizing the specific lines of an image by binarization and uses all the pixel points in an image to find the parameters of a straight line. Each point obtains its parameter value by one-to-many projection to the parameter space, and parameter values are generated by all the points to find the most obvious line parameter in the space. The Generalized Hough Transformation does not have much influence on image distortion or pixel loss and is not limited to the recognition of straight lines. In fact, it can be used to recognize any shape. However, it cannot make an efficient classification, and the outcome is not stable enough for computer vision. The parameter space used to search for the target can become very large in a short time, and a lot   of storage space is needed. However, the Hough transformation can solve the above-mentioned problems, and it was chosen for this study.
It can be seen in Table 2 that after improvement, the Hough transformation method can quite accurately locate the mass point of the workpiece to be gripped. The pixel error for the mass center was also reduced by 3.09 pixels. However, the catch rate, when compared to the method before improvement, was reduced by 26.27%. However, the time for image acquisition (before the robotic arm grips the workpiece) is more than adequate for the Hough transformation to determine the mass center of the workpiece to be gripped. Therefore, image recognition, as performed by our improved Hough transformation method, was used in this study.

Conclusion
In this study, communication between the "image" and the robotic arm was proposed to reduce preexisting errors in location. Integration of the robotic arm and image identification allowed automatic judgment to be made without the need for an extra manipulator. Correction of the inherent errors and accurate calibration was achieved via feedback vision. This allowed the endpoint pose position of the arm to be fine-tuned, and this significantly improved its accuracy. Measurements were made with the highly stable Hough transformation  method rather than by the in-depth mechanical learning methods used before. Hough transformation works without the need to build a huge database of training data and saves a lot of time. Results showed that changing the parameters to increase the Hough transformation accuracy not only reduced image pixel error but also increased the precision with which the mass point position could be determined and gave better gripping performance.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was supported by the Ministry of Science and Technology under Project number: MOST 107-2218-E-167-001.


extfieldsofstudy:
Computer Science
Engineering

text:
Wearable Travel Aid for Environment Perception and Navigation of Visually Impaired People

This paper presents a wearable assistive device with the shape of a pair of eyeglasses that allows visually impaired people to navigate safely and quickly in unfamiliar environment, as well as perceive the complicated environment to automatically make decisions on the direction to move. The device uses a consumer Red, Green, Blue and Depth (RGB-D) camera and an Inertial Measurement Unit (IMU) to detect obstacles. As the device leverages the ground height continuity among adjacent image frames, it is able to segment the ground from obstacles accurately and rapidly. Based on the detected ground, the optimal walkable direction is computed and the user is then informed via converted beep sound. Moreover, by utilizing deep learning techniques, the device can semantically categorize the detected obstacles to improve the users' perception of surroundings. It combines a Convolutional Neural Network (CNN) deployed on a smartphone with a depth-image-based object detection to decide what the object type is and where the object is located, and then notifies the user of such information via speech. We evaluated the device's performance with different experiments in which 20 visually impaired people were asked to wear the device and move in an office, and found that they were able to avoid obstacle collisions and find the way in complicated scenarios.

I. INTRODUCTION
A CCORDING to the World Health Organization (WHO), there are about 285 million visually impaired people all over the world [1]. Visual impairment makes it challenging for them to perform autonomous navigation and environment perception in unfamiliar environments [2], [3]. Visually impaired people usually asked white canes or guide dogs for help when detecting obstacles over the past years. However, they may still suffer injuries from hanging obstacles, such as scaffoldings and portable ladders, as white canes and guide dogs can only detect obstacles at heights up to their chests [4], [5]. Recently, Electronic Travel Aids (ETAs) [6] utilizing advanced sensing techniques have greatly improved the travelling experience of visually impaired people. However, the ultrasonic sensor based ETAs are poor at obstacle identification due to the wide beam angle of ultrasonic sensors, and the laser sensors based ETAs are expensive, heavy and have high power consumption, which make them unsuitable for wearable applications [4], [7]. Although vision based ETAs (e.g. monocular camera J. Bai   based [8], stereo camera based [9], [10] and Red, Green, Blue and Depth (RGB-D) camera based [11], [12]) have been widely used for assisting visually impaired people to avoid obstacles and find obstacle-free paths, some problems still exist. For example, a chair may be considered as an obstacle when a blind person is looking for a seat, or no path will be found to enter a room with closed doors. Whereas suppose the object recognition technique be adopted, the visually impaired people can find the chair to sit or open the door to enter the room. Therefore, there is a need for developing new assistive devices to help visually impaired people in navigation and environment perception. This paper proposes a wearable assistive device (see Fig. 1) for visually impaired people's navigation and environment perception. The presented device detects the ground accurately and rapidly by leveraging an adaptive ground height segmentation algorithm and utilizing the ground height continuity among adjacent frames. With the detected ground, the optimal walkable direction can be computed and corresponding beep sound will be played out for blind navigation. Meanwhile, we adopt a lightweight Convolutional Neural Network (CNN) to semantically categorize the detected objects in the RGB image and then extract the object contour in the depth image to identify where the object is located. All the semantic information is finally converted to speech to inform the user about the perceived surroundings.

II. RELATED WORK
This research builds on important related work of ground detection and object recognition. Because the vision based assistive devices have advantages over the ultrasonic or laser sensor based devices as mentioned above, this section focuses the relevant works on vision based ground detection and object recognition. General ground detection algorithms can be divided into scene segmentation based algorithms and surface normal vector estimation based algorithms.
Segmentation: Saitoh et al. [10] proposed a mean-shift algorithm to fit the ground plane. If the angle between the fitting ground and the horizontal plane is less than a threshold, the ground will be taken as a traversable area. However, the algorithm is too sensitive to the threshold. Rodriguez et al. [14] proposed a Random Sample Consensus (RANSAC) [15] based ground plane detection method and the potential obstacle was represented by polar grid. The ground detection error of this method reaches more than ten percent, leading to high obstacle detection errors. The wearability of the system and the integration of its main components into smaller devices should also be improved.
Surface normal vector estimation: Koester et al. [16] presented a gradient and surface normal vector based detection algorithm to compute the accessible sections effectively even in crowded scenes. However, the success rate of detection heavily relies on the quality of 3-Dimensional(D) reconstruction process. Bellone et al. [17] estimated the normal vectors to a local surface via Principal Component Analysis (PCA) and generated an unevenness point descriptor to detect traversable and non-traversable regions. The performance of the method was greatly affected by the search radius, which limits its applications in practice. Aladren et al. [7] used depth and color images to detect the ground in longer distance. Although their system achieved high ground detection accuracy, it is too computational expensive (2 frames per second) for visually impaired people to navigate in real time. Similar to our work, Imai et al. [18] proposed a ground detection approach which considered both ground height and normal vectors. However, since the ground height was computed only by one column data in depth image, it is prone to error if an obstacle happened to exist in this column. In contrast, the ground height in our approach is computed more robustly and rapidly through weighting the ground heights in previous and current frames.

B. Object Recognition
Tapu et al. [21] presented a real-time obstacle detection and classification system for safe navigation of visually impaired people. The system used Scale Invariant Feature Transform (SIFT) and Features from Accelerated Segment Test (FAST) features to extract points of interest, and used Support Vector Machine (SVM) and Bag of Visual Words (BoVW) to classify these points. Although the system achieved 90% precision of classification, the information about object distance is not available, and this will result in misguidance to blind individuals. Lee et al. [11] proposed a robust depth-based obstacle detection system to obtain obstacle information which contains the distance, while it did not employ object recognition techniques, leading to poor environment perception functions.
Since AlexNet [22] won the ImageNet Challenge: ILSVRC 2012 [23], CNN based object detection methods have become unprecedentedly popular. Although these methods have higher detection accuracy, they usually have large network sizes and high hardware performance requirements. Kaur et al. [24] proposed a faster Recursive Convolutional Neural Network (RCNN) based scene perception system for visually impaired individuals. The system is able to provide the obstacle category and distance information. However, the distance information is obtained only through a single-line laser, making it unavailable or prone to error for some obstacles. Tapu et al. [25] proposed a DEEP-SEE framework that used both computer vision algorithms and CNN to detect objects encountered during navigation. Although its recognition accuracy is satisfactory, the heavy computation load makes it difficult to be implemented on a smartphone.
Lightweight CNN based 2-D object detection algorithms [26] [27] [28] [29] [30] have made great progress in recent years. They usually provide object category and location information in 2-D images, while there is still a lack of distance information. As a result, an object painted on the ground may be considered as an obstacle, and the visually impaired people will be misled.
To overcome the above limitations, we propose a 2.5-D object detection method that can provide the object category, distance and orientation information to make visually impaired people travel more easily.

III. SYSTEM DESIGN
As shown in Fig. 2, the system first acquires RGB image and depth image from an RGB-D camera, and also obtains the camera attitude angle from an Inertial Measurement Unit (IMU) attached on the camera. Then based on the above sensor data, a time-dependent adaptive ground detection algorithm is performed to detect ground. Next, an optimal walkable direction search algorithm is employed to find the direction that visually impaired people can follow. This direction will then be converted to beep sound to give the users instructions. When the user wants to know the surroundings, he can double tap the smartphone screen to trigger the 2.5-D object detection function. The surrounding information including category, distance and orientation of all obstacles will be provided to the user through speech. The algorithms will be described in the following sections in detail.

A. Time-Dependent Adaptive Ground Detection
Pointcloud reconstruction: As shown in Fig. 3, the camera coordinate system X c Y c Z c is centered at the camera, and the positive Z c −axis, Y c −axis, and X c −axis are defined as the camera's facing direction, up direction and left direction respectively. The world coordinate system X w Y w Z w is centered at the camera coordinate system center, the positive Z w −axis, Y w −axis, and X w −axis are the user's facing direction, vertically upward direction, and left direction respectively. Both coordinate systems are the left-handed Cartesian coordinate systems. The pixel value of point p(u, v) in the depth image represents the distance between point P (x, y, z) and the camera, which is equal to z. With the camera attitude angle measured by the IMU, the corresponding 3-D pointcloud in the world coordinate system can be calculated through: where K is the camera intrinsic parameter matrix, point (x w , y w , z w ) is the reconstructed point in the world coordinate system, α is the camera pitch angle, γ is the camera roll angle. Fig. 3. Coordinate system transformation.
Coarse ground fitting: As shown in Fig. 4, the initial ground height threshold T Y roi is calculated adaptively using the OTSU algorithm [31] in the current frame. Since the change of ground height in two adjacent frames is usually limited, the ground height T Y pre of the previous frame is used to reduce the perturbation of other planes (e.g. desk, sofa) (see Fig. 9). The final ground height threshold is computed as: where λ, µ are the weights.
Due to the inherent limitation of the depth camera, the depth accuracy always drops down with the increase of distance. Besides, the obstacles that are too far away from the person do not need to be considered. Therefore, only the points within a threshold T Z are used in order to reduce computation cost. By making use of the ground height T Y and the distance threshold T Z, the 3-D points for fitting the coarse ground can be computed as: Then the coarse ground is fitted with RANSAC algorithm [15], and represented as: Ground refinement: The normal vector n ground of coarse ground plane can be obtained directly by Eq. 4, and the ground pitch angle φ will then be computed through: where n xoz is the normal vector of plane X w OZ w . According to the ground pitch angle and the empirical slope angle, the coarse ground will be classified as one of the four types: horizontal, upslope, downslope and non-ground. If it is non-ground, then the visually impaired people will be directly informed that they cannot move on; otherwise, the coarse ground will be refined with the unevenness tolerance σ through: where dist(p, F init ) is the distance from point p to the coarse ground, F is the final 3-D point cloud of refined ground.
Finally, the refined ground height H is obtained through: and it will be used in next frame.

B. Optimal Walkable Direction Search
If no ground is detected, the system will directly inform the visually impaired people and stop proceeding the optimal walkable direction search algorithm. Otherwise, the optimal walkable direction is calculated as follows.
Since the walkable direction relies on the detected ground, only the 3-D points within the ground plane need to be considered for walkable direction search. These points are selected by: where F is the 3-D points on ground plane, H is the ground height, A, B, C, D are the plane parameters defined in Eq. 4, and ε is a constant for preventing the person from being collided with overhanging obstacles.
Then the 3-D points P are projected onto the plane X w OZ w (see Fig. 5). The nearest points in all sectors (each sector is represented as the sub-region in Fig. 5 and the angle of each sub-region is 0.5 • ) can be easily obtained. The award of each sector is computed as: where θ is the angle of a sector, w sw is the passable width (greater than the person's body width), z i is the z−axis coordinate values of the nearest point in sector i, α and β are the weights, and N is the total number of sectors. This award function ensures that the person moves toward the direction with smaller turn angle and longer traversable distance.
Next, the optimal walkable direction is obtained by: where τ is a distance threshold, z arg max i (award[i]) is the nearest distance of the sector with the maximum award, i max is the index of the sector with the maximum award, and N is the total number of sectors. If z arg max i (award[i]) is less than a small value τ , it is considered to have a very large risk of collision with obstacles. In that case, the optimal walkable direction does not exist, and the system will inform the users to turn left or right with a large angle to search a walkable direction. If z arg max i (award[i]) is larger than τ , the optimal walkable direction is that corresponding to the sector with the maximum award. If the turning angle is very small (e.g. |γ| ≤ 5 • ), the system will directly inform the users to go straight to prevent them from vacillating to the left and right.

C. 2.5-D Object Detection
CNN based object detection: As the MobileNet V2 [32] achieves relatively better results than other networks [26] [28] [29] [30] by using depth-wise separable convolutions, and has lower computational complexity and smaller model size, it is utilized here for object detection. The training is implemented on the COCO dataset, which includes 91 classes, such as person, car, bus, chair, etc. These classes are enough for helping visually impaired people have a general perception of surroundings. However, 2-D object detection cannot be directly used for blind navigation due to the lack of object distance information. For example, if an object painted on ground is recognized as an obstacle, it will lead visually impaired people to make an incorrect decision. Therefore, we also use the depth image based object detection to solve the above problem.
Depth image based object detection: The detected ground (see Section III-A) is firstly removed from the depth image. Then the close morphological processing 1 is performed to merge the small objects. Next, the external contours of the obstacles are extracted and their areas are computed. If the area is less than the threshold S, the corresponding obstacle will be merged into its nearest obstacle or taken as noise; otherwise, the obstacle location can be obtained as: 1) Compute the moment of each contour, and the centroid of the contour can be obtained according to the zeroorder and the first-order moment; 2) According to the contour centroid and the camera intrinsic parameter matrix K in Eq. 1, the obstacle location is represented as (θ pitch , θ yaw , z): where center(x, y) is the contour centroid, and z center is the depth value of contour centroid. Combination: The objects obtained by MobileNet V2 in the RGB image can be easily mapped to the depth image according to camera calibration matrix [R, t]. Then the intersection area C between the mapped area A and the detected contour B can be calculated (see Fig. 6): If S C max(S A ,S B ) (S i represents the area of region i) is greater than a threshold ζ (e.g. 0.7), it means the mapped area and the detected contour is the same object. Then the object distance is the minimum non-zero depth value within the intersection area C, and the object orientation relative to the user can be obtained with Eq. 11. This way, the key surrounding information including obstacle category, distance and orientation can be provided.

D. Audio Feedback
Navigation feedback: The beep sound is used to provide the navigation information to the user. When the visually impaired person encounters an obstacle that blocks his way, the system will keep beeping to alarm him that he cannot go straight. In that case, he should turn left or right to search the optimal walkable direction. When the beep sound stops, the user can continue moving.
2.5-D object detection feedback: When the visually impaired person walks in a relatively complicated environment, the 2.5-D object detection function can be activated by double tapping the smartphone screen, and the results will be converted to speech via a text-to-speech module and broadcast to the user. An example of 2.5-D object detection feedback is shown in Fig. 7, the object category (e.g. chair), distance (e.g.

IV. EXPERIMENTAL RESULTS AND DISCUSSIONS
The performance of the proposed system is evaluated in real scenarios. Firstly, the ground detection is tested since it plays an important role in the whole system. Then, we recruit 20 visually impaired persons to test the navigation and object detection function. We follow the protocol approved by the Bejing Fangshan District Disabled Persons' Federation for recruitment and experiments. All the people who participated in this experiment approved that the results (including data, images and videos) can be published with anonymity. Finally, the computational cost is evaluated to test the real-time performance.

A. Experiments on Ground Detection
To evaluate the proposed ground detection method quantitatively, we manually labelled a random sample of 1000 images captured in indoor scenario to get the ground truth. Fig. 8 shows the precision of the detected ground with different distance and different Intersection over Union (IOU) percentages. The IOU percentage is computed as: where N ∩ is the number of overlapping pixels between the detected ground and the ground truth, N detected is the number of detected ground pixels, and N groundtruth is the number of ground truth pixels.  As shown in Fig. 8, comparing with the RANSAC algorithm, the proposed method has a higher precision in all scenarios. Some comparative examples are shown in Fig. 9. By taking the advantage of the ground height continuity among adjacent frames, the proposed method is more robust to the interferences of other planes (such as sofa, table) than the RANSAC algorithm. It provides a solid base for the following object detection.

B. Experiments on Real Scenarios
Navigation task: Three paths in an office (see Fig. 10) are selected to evaluate the navigation performance. There are many obstacles including static and moving obstacles in these paths. Then 20 visually impaired people who were not familiar with these environments are asked to navigate following these paths with the help of either a white cane or the proposed device. They are trained for about 10 minutes to get familiar with our system. The average walking time and number of collisions with obstacles are recorded (see TABLE I). The visually impaired people using our proposed system spend less walking time than using the white cane. This proves that the proposed system has a higher navigation efficiency in unfamiliar environments.  As shown in TABLE I, the users have more collisions with obstacles when using a white cane. This is because they usually use the white cane to detect the obstacles on the ground instead of those hanging in the mid-air, and it is observed that almost all studied users collide with the desk or the hanging objects on the path C→D and E→F (see Fig. 11). Whereas with the proposed device, they are able to avoid collisions with those obstalces. This proved that the proposed system is more secure for visually impaired people.
2.5-D object detection task: Some complicated scenarios are designed to test if the 2.5-D object detection can help visually impaired people to travel more efficiently. Such examples include the scenario where the way is blocked by two chairs (see Fig. 12) and crowded with other obstacles. With only the navigation function in our proposed system, the user cannot Collide with the desk Collide with the hanging object Fig. 11. Collisions when the visually impaired people use a white cane.
find the moving direction after multiple searches and will turn back. However, if he activates the 2.5-D object detection and perceives more information about surroundings, he is able to move the chair and search the moving direction. This prevents the visually impaired people from taking more detours and improves their environmental perception abilities. This shows that the 2.5-D object detection indeed helps visually impaired people to travel more efficiently and brings a better travel experience to them.

C. Computational Cost
All algorithms are implemented on a smartphone with Qualcomm Snapdragon 820 CPU 2.0 GHz and RAM of 4 GB. The average computational time of the proposed system is calculated and shown in TABLE II. The images acquisition, ground detection, optimal walkable direction search and 2.5-D object detection cost about 0.66 ms, 13.53 ms, 7.19 ms and 114.13 ms respectively. The total time that all algorithms cost excluding the 2.5-D object detection is about 27.17 ms on average. Since the 2.5-D object detection is only activated when the visually impaired people want to know the surrounding information or enter a complicated scenario, the proposed system is able to provide real-time assistance for users' daily traveling. V. CONCLUSION This paper presents a wearable device which provides navigation and object detection assistance for visually impaired people. To provide reliable navigation, a ground detection algorithm that uses the ground height continuity between two adjacent frames is presented. Then an optimal walkable direction search method is developed to determine the moving direction. To improve the environmental perception ability of visually impaired people, a 2.5-D object detection function is presented. Audio feedback is used to inform the visually impaired people for both the navigation instructions and object detection results. Experimental results show that the proposed system can help visually impaired people travel efficiently and brings a better traveling experience for them. In future, the small-size obstacle detection will be considered.