Proving ground, or on-track testing has been an essential part of testing and validation process for connected and autonomous vehicles (CAV). Several world-class CAV proving grounds, such as Mcity at the University of Michigan and The Castle of Waymo, have already been built, and many more are currently under construction. In this paper, we propose the first optimization approach to CAV proving ground designing and refer to any such CAV-centric design problem as “Xcity” to emphasize the enormous investment, the multi-dimensional spatial consideration, and the immense construction effort emerging globally. Inspired by the recent progress on traffic encounter clustering, we further define “road assets” as fundamental building blocks and formulate the whole design process into nonlinear optimization problems. We have shown that such framework can be utilized to adaptively generate CAV proving ground designs with optimized capability and flexibility and can further be extended to evaluate an existing “Xcity” design.
Keywords: connected and autonomous vehicle, testing and validation, proving ground design, traffic encounter, optimization
Development of self-driving technologies has drawn significant public attention with major concerns on safety issues. Extensive efforts have been made on testing and validation (T&V) of both autonomy and safety of connected and autonomous vehicle (CAV) systems before deployment. The "V-model" for CAV systems [IgorP2015Vmodel] illustrates the T&V process from component level to traffic level. An approach for such process, visualized as a T&V pyramid [zszalay2016pyramid], has identified proving ground testing as an essential level required by CAV systems. CAV proving ground is a reserved area designed and constructed to simulate real-world traffic environment. The fully controllable environment of proving grounds provides CAV testing and validation with foreseeable risks and flexible fidelities. Those advantages have led to the establishment of several world-class CAV proving grounds: Mcity at the University of Michigan [mcity], "The Castle" of Waymo [castle], Uber’s ALMONO [almono], Smart Mobility Advanced Research and Test Center (SMART) by Transportation Research Center (TRC) [trcsmart], Security Smart Mobility Analysis and Research Test (SMART) Range in Israel by HARMAN [harman], American Center for Mobility (ACM)’s test facility at Willow Run in Michigan [acm], etc. Since massive time and resources have been invested in CAV proving ground construction, it is natural and essential to aim at effective and efficient designs. However, design theory of CAV proving grounds is rarely seen in the existing literature.
When designing an efficient and effective CAV proving ground, one should make clear the testing capability required by CAV T&V starting from elementary system integration, such as path following, till high-fidelity traffic simulation, such as CAV tests. Due to the infeasiblity of complete testing on CAV systems [philipK2016challenge], it is natural for one to break multi-vehicle-multi-event and long-term CAV tests into a variety of small and distinguished use cases which represent typical usage scenarios for autonomous driving [Wachenfeld2016]. Such typical use cases or scenarios are similarly identified in [PATH2014, NHTSA] as behavioral competency, the ability of a CAV to operate in traffic conditions which it will regularly encounter; these conditions include keeping the vehicle in the lane, obeying traffic laws, following reasonable etiquette, and responding to other vehicles, road users, or commonly encountered hazards. Due to the highly complicated and stochastic nature of real traffic, it would be impractical to require tests for all possible use cases. Therefore, it is reasonable for one to identify minimum behavioral competencies as required. Such efforts lead to CAV capability checklists which will be practical for one to use as evaluation criterion for CAV systems [PATH2014, TORC]. Thus, a CAV proving ground is considered both practical and effective when it provides wide coverage of identified behavioral competency tests within constrained space. To the best of our knowledge, the literature provides no such systematic approach that employs optimization models to maximize the testing capability of a CAV proving design. This work contributes to the formulation of a systematic proving ground design that cognizantly addresses the CAV evaluation challenges to assist the design of CAV proving grounds in terms of maximizing testing capability and flexibility.
Challenges of systematically designing an effective and efficient CAV proving ground come from two aspects: a) it is not clear how to classify and extract typical driving scenarios from a large scale naturalistic driving data; b) it is not clear how to map those scenarios into space-constrained proving grounds. To address these problems, we propose the development of Xcity: a mathematical approach to design CAV proving ground road layout based on traffic encounter theory[trafficprimitive, wenshuo, sisiLi] and nonlinear optimization methods. We applied and illustrated our framework with a minimal example. Then we extracted driving scenarios from Mcity road map and from clustered driving scenarios by [wenshuo], thereby exploring the scalability of our approach.
3 Xcity Design Framework
The underlying idea of Xcity design framework is to generate a CAV proving ground road map which provides maximum coverage of typical CAV driving scenarios within constrained space and maximum flexibility of scenario transitions. In this section, we describe the systematic approach in terms of scenario extraction and representation, scenario selection, and placement optimization.
As mentioned in introduction section, we evaluate the capability of a CAV proving ground based on the CAV driving scenarios it supports. It is worth mentioning that this paper focuses on scenarios involving at least two vehicles. The reason behind this is that single-vehicle tests such as lane following, traffic rule obeying, and intersection navigations can be triggered by simply removing all the other road users from the corresponding multi-vehicle scenarios, e.g. car-following including Stop & Go & Emergency Stop and 2-way stops obeying with cross traffic [TORC, PATH2014]. Thus, we argue that the evaluation requirements of single-vehicle tests can be trivially incorporated whenever more complicated multi-vehicle scenarios are conductable.
Therefore, it is essential to first develop a systematic way to identify, extract, and represent typical multi-vehicle CAV usage cases from a large-scale driving data. Such task is non-trivial since the real-world traffic is a vast, stochastic, and dynamic cyber-physical system consisting of large number of road users, non-road users, and stochastic environment factors [trafficprimitive]. The method developed in [sisiLi, wenshuo, trafficprimitive] allows the processing of massive multidimensional traffic data to extract traffic primitives and traffic encounters. In this context, traffic primitives, which are the most informative traffic features from the dataset, are then used as principal compositions of the entire traffic [trafficprimitive]. Traffic encounters, which represent the scenario where two or multiple vehicles are spatially close to and interact with each other on roadways [wenshuo], allow higher fidelity representation of sophisticated driving scenarios.
Our design approach aims to map such scenarios into the given space in a way that maximizes CAV evaluation capability. To enable tests on a particular scenario, we map the group of roads and intersections which supports the realization of traffic encounters into the proving ground design. Additional factors such as other road users, traffic control devices, and control commands are then embedded on demand according to the test specifications. For example, one could add two stops signs and deploy one other vehicle at an X-intersection to support the scenario of 2-way stop-obeying with cross traffic. With the same physical layout, one could also evaluate any scenarios of the form left-turn-on-green with incoming traffic. Therefore, we consider the combination of roads and intersections as basic construction element in our design formulation and refer to such combination as road assets hereafter. Once selected and mapped into the design, a road asset will be treated as a 2D rigid body afterwards and all corresponding testing scenarios it supports (based on traffic-primitives analysis) will be enabled. Furthermore, we introduce a value metric to evaluate road assets since the number of supported scenarios for different assets could vary. For a road asset , its value comes from both its versatility to support multiple scenarios and the universality of the supported scenarios.
Given a preprocessed set of road assets , value measure , and constrained construction space , we propose the following constraints, definitions, and objectives to model an Xcity design problem.
Roads in any single asset cannot overlap with roads in other assets.
Road assets should completely reside within the constrained space .
A subset, , of is considered feasible if all road assets in can be mapped into simultaneously under the above constraints.
Two road assets and are considered directly connectible if it is possible to construct a straight transition road to connect them without colliding with any other assets. The number of directly connectible asset pairs in a road asset set is referred to as direct connectivity hereafter and denoted as .
The feasible subset with highest total asset value is selected and mapped in .
The assets in are arranged within such that is maximized.
Constraint 1 guarantees that no test fidelity is lost due to the physical design. Each scenario test should be conducted in its originally dedicated road asset with unaltered roads and intersections. Merging distinguished assets would generate undesired and inaccurate physical environment, which attenuates the fidelity of the simulated scenario. It would consequently cause non-measurable degeneration to the confidence of public street deployment gained from on-track tests.
As previously mentioned, one criteria which we use to evaluate a CAV proving ground is its coverage of behavioral competency (BC) tests. BC is defined as a list of independent scenarios that are expected to be handled safely by a sophisticated autonomous driving system. Since BC implies no connections between different scenarios, we treat the mapped road assets independently and simply accumulate their value measures to form Objective 1. On the other hand, Proving ground test and validation, a comprehensive and multi-integration-level process, calls for not only a wide range of independent scenario-based tests but also cascaded and potentially randomly ordered scenarios to assess CAV systems’ robustness and performance on enduring operations. Thus, the success in emulating scenario sequences further boosts the flexibility and fidelity of the proving ground. Tests involving multiple independent but connected scenarios highly demand that the CAV proving ground be flexible on navigating among different scenarios. In our case, such flexibility will be ensured by transition roads built to connect different road assets. Although roads with arbitrary shapes and lengths would highly likely to serve as transitions between any pair of road assets, it is worth looking for a design with more directly connectible asset pairs. Therefore, we only consider straight transitions for evaluation purpose. The flexibility and fidelity measure depending on possible straight transition roads are defined as Objective 2. Note that the the outcome of this paper will be a high-fidelity graph representation of the road map structure. Detailed information on road construction, such as exact road curvatures, lane widths, and sidewalks, needs to be added to our design, so as to achieve a final construction blueprint. Further discussion on blueprint rendering is out of scope and is omitted.
In this section, we describe the mathematical formulation and optimization models of our approach. With road assets pre-produced and represented, the design process is decomposed into two phases. In phase 1, the feasible road asset subset with highest total asset value is selected from all pre-produced road assets . In phase 2, the placement of the pre-selected within given space will be optimized in terms of cross-scenario transition flexibility. Both two phases are modeled into nonlinear optimization problems.
4.1 Road Asset Representation
Traffic encounters are clustered as GPS location sequences of involved vehicles. Thus, given a traffic encounter, one natural and direct approach to retrieve the essential road asset is to directly locate the traffic encounter on the world map and extract surrounding roads and intersections. To map the extracted road assets into the given space, we will need to define a geometric descriptor for their structure and size. On one hand, such descriptor should be rigorous enough for accurate localization and space-reserving for the corresponding asset. On the other hand, the descriptor should be concise enough to truncate unnecessary computations during optimization. Two of the most recognized map representations for autonomous driving are Lanelets [lanelets] and OpenDRIVE [opendrive]. Both of them are efficient and detailed enough for semantic environment modeling and motion planning. However, our approach only aims at the placement and transitions of road assets regardless of their semantic meanings or navigational information. These knowledge is already contained in the supported scenarios and can be applied on demand after the design is finished. Therefore, both Lanelets and OpenDRIVE contain more information than needed and create unnecessary computations.
Considering both representativeness and simplicity, we use graphs to model the geometric structure of road assets (see examples in Figure 1). Such descriptor is inspired by the map format from OpenStreetMap (OSM) [OpenStreetMap], a collaborative project to create a free editable map of the world [wiki:osm]. We mark the intersections, road ends, and way-points along curved roads with nodes . Then, we connect each node pair , if necessary, with a line segment to model one road segment in the real world. We refer to those line segments as the internal segments of road assets. We ignore nodes that are unnecessary or trivial, thereby maintaining the whole structure to introduce the least possible internal segments and save computation efforts. Typical trivial nodes include nodes at X-intersection centers (see Figure 1.e). Notably, the number of internal segments turns out to be one major factor contributing to model complexity.
4.2 Phase 1: Most-valuable Feasible Subset Selection Model
4.2.1 Single-Asset Constraint Set (Sacs)
We map road assets into the given space by mapping their corresponding graph representations while retaining the relative placement of the nodes. This is done through the introduction of a set of 2D geometric constraints on the nodes of which the original node placement is the only realization. A non-trivial start point would be mapping a triangle defined by the original vertex positions into space at without deformation or flipping. We propose the following constraint set , referred to as triplet constraint set, which retains the triangle’s original shape and size. Note that the superscript number refers to the node index within a single asset. The subscript ’o’ refers to the original location of nodes on world map.
, the orientation test, is positive if are placed counter-clock-wise (CCW), or negative if clock-wise (CW), or zero if collinear (LNR) [orient]. See Equation 7 for its definition. The variables in constraint 1-3 are introduced as slack variables for equality constraints to avoid numerical issues in solvers. The sum of all squared slack variables will be include in the objective to be minimized.
The relative placement of all nodes in road asset can be retained by applying repeatedly to nodes in . We refer to the constraint set which regulates the relative placements of all nodes within a single asset as Single-Asset Constraint Set (SACS), denoted by (see constraint set 8).
4.2.2 Cross-Asset Constraint Set (Cacs)
Another requirement of road asset mapping is that internal segments of different road assets should not collide or intersect. To formulate such requirement into constraints, we first define the intersection test for straight line passing points and line segment formed by points as Equation 9.
is then positive if doesn’t intersect with and non-positive otherwise. See Figure 3 for illustrative examples.
Given two internal segments and , we denote the straight lines passing and as and respectively. Then doesn’t intersect with () if Constraint 10 holds.
To simplify the constraints, we skip considering the corner case where and are collinear but not overlapped. Such case can be trivially made feasible under our constraints by rotating or moving one of the line segments by a tiny amount. It is worth mentioning that, since strict inequalities are not permitted in most solvers, we implement constraint 10 using non-strict equalities with a constant offset
. Auxiliary binary variablesare introduced as well to decompose the logic OR statement into a set of AND statements for each instance of constraint 10. Constraint 10 is finally formed as constraint set 11. We refer to such constraint set that eliminates the collision of and as Cross-Asset Constraint Set (CACS), denoted as .
4.2.3 Phase 1 Optimization Model
With both SACS and CACS defined, we can formulate phase 1 in a single optimization model (see Expression 4.2.3). For each subset , let be the collection of all mapped nodes in space from . Let
be a column vector containing all slack variables introduced in distance equality constraints. Letbe the set of binary variables introduced in CACSs. |l| A’⊆A,X,b_CACSδ^Tδ- υ(A’) K_sa(a_τ) τ=1,2,…,NA’ K_ca(ℓ_p,ℓ_q) ∀ ℓ_p∈a’_i, ℓ_q∈a’_j⇒a’_i, a’_j∈A^′, i≠j Notably, we factor out the selection of asset subset during implementation and manually assign the road assets to be mapped. When global minimum of the remaining objective is reached, all distance Constraints 1-3 are strictly satisfied. The corresponding will be a feasible mapping of into .
4.3 Phase 2: Transition Flexibility Optimization Model
With the most valuable feasible subset selected in phase 1, we will find a placement of all containing road assets such that the direct connectivity, , is maximized.
4.3.1 Asset Transition
In the graph representation of a road asset, we define nodes that are owned by only one internal segment as boundary nodes (for instance, , , in Figure 1.a) and the rest as internal nodes. To avoid causing bias to the road layout and achieve high-fidelity scenarios, we only allow road assets to be connected to each other using boundary nodes and skip the internal nodes when considering transition road construction. Two road assets are considered directly connectible if it is possible to construct at least one straight transition road between them using only boundary nodes without intersecting any internal segments of any road asset. See Figure 4 for a mini example.
4.3.2 Phase 2 Optimization Model
Given a the pre-selected best feasible asset set from phase 1, the objective is to maximize subject to both SACS and CACS. Given a road asset with nodes , we propose the following definitions. Note that the superscript * indicates that the variable is pre-selected in phase 1 and kept constant in phase 2.
: set of binary decision variables indicating the direct connecting status of and .
: set of boundary node in .
: set of transition road candidates .
: binary variable which indicates the selection status of .
: set of internal segments that don’t have common node with by construction.
With SACS and CACS defined in Equation 8 and 11, we can write the complete optimization model as Expression 4.3.2. |l| X, b_dc, b_tr, b_CACS δ^Tδ- sum(b_dc) ∣ A^* K_sa(a_τ) τ=1,…,NA^* K_ca(ℓ_p,ℓ_q) ∀ ℓ_p∈a^*_i, ℓ_q∈a^*_j⇒a^*_i, a^*_j∈A^*, i≠j ∑_p,qb_tr(x_i^(p), x_j^(q))≥1 ∀a^*_i, a^*_j∈A^*⇒i≠j, b_dc(a^*_i, a^*_j)_true K_ca(ℓ_tr,ℓ_χ) ∀ℓ_tr⇒b_tr(ℓ_tr)_true, ∀ℓ_χ∈L_χ(ℓ_tr) Constraints 4.3.2 and 4.3.2 ensure that the feasibility of the mapping as identically done in phase 1. Constraint 4.3.2 ensures that for each pair of assets to be connected, at least one transition road candidate connecting their boundary nodes should be selected. Constraint 4.3.2 ensures that each selected transition road candidate should not intersect with any other internal segment unless they share a node by construction. In Figure 4, to is one example of such exception.
5 Optimization Approaches
Based on our formulation, the resulting model is NP hard, essentially due to the extensive use of binaries to ensure the validity of the asset shapes and their combined designs. From (4.2.3) and (4.3.2), it is clear that even though the final objective function in (4.3.2) is arguably convex, some of the constraints in both (4.2.3) and (4.3.2) are not. As such, the model is of mixed integer programming (IP) type and the natural approaches to numerically deal with it take on some form of branch-and-bound (BNB) methods. Loosely speaking, the computation procedure of BNB requires lower- and upper-bound governing how variable branches are determined to obtain integer solution at each numerical iteration [lawler1966branch]. Modern solvers employ various numerical scheme deemed best to solve these branch-level problems, e.g. BMIBNB use linear relaxations to solve for the bounds from the branch problems [lofberg2004yalmip]. Furthermore, other solvers have gained improved efficiencies from more sophisticated techniques to assist the branching, e.g. SCIP [scip1] employs branch-cut-and-price (BCP) while BARON [baron1] uses branch-and-reduce (BAR). In what follows, we describe the performance of these three solvers in regard to our formulation.
BMIBNB is a YALMIP-based BNB technique that can be easily implemented using standard computing environment Matlab [lofberg2004yalmip]. YALMIP itself is an interface program that is originally designed to allow Matlab to solve semi definite programming problems using external solvers including SeDuMi and Mosek, but its implementation has currently included wider array of choices including fmincon, SDPT3, Gurobi, and CPLEX which can be tuned easily [lofberg2004yalmip] to obtain the numerical bounds. In dealing with nonlinear IP, BMIBNB uses linear relaxation and convex envelope approximations. As such, it trades the accuracy of solution optimality for the efficiency of computation. With this view, one shall notice that when the original problem is non-convex and high-dimensional, then the performance of BMIBNB suffers from yet-to-converge bounds. Unfortunately, our Xcity design problem is one such program and thus, as expected, the corresponding numerical experiment advises to using more customized solvers for our formulation.
SCIP is a BNB-based solver for integer programming enhanced with BCP framework which was designed to handle large-scale discrete optimization problem. The cuts and variable can be generated dynamically through the branching process to explore the direction of optimality. This dynamic ‘exploration’ capability is achieved through the implementation of column generation technique, which further implies SCIP has adopted the large-scale Dantzig-Wolfe decomposition method [scip1, scip2]. Furthermore, the most recent version has included Bender’s decomposition, allowing the solver to efficiently deal with problems with block-angular characteristics [geoffrion1972generalized]. Our preliminary numerical experiment shows that the SCIP’s inherent column generation does not yield an appealing performance due to the insurmountable growth of the constraint size in both the number and shape-complexity of the assets.
BARON (Branch-And-Reduce Optimization Navigator) is designed to solve a generalized nonlinear programs including mixed integer nonlinear problem [baron1, baron2]. The underlying idea of this technique is that after branching, a numerical analysis is performed to reduce the potential optimal regions which in turn allows a faster computation, and if conducted appropriately, global optimality assurance [kilincc2018exploiting]
. This notion of ‘appropriateness’ is achieved by adapting the cutting plane and probing techniques from mixed integer linear programming to mixed integernonlinear problems. We refer interested readers to [kilincc2018exploiting] for extensive computational experiment investigating the effectiveness of BARON scheme to achieve optimality. Our numerical results show that BARON indeed yields to much better performance both in terms of improved computational efficiency and solution quality compared to BMIBNB and SCIP. As such, we will use BARON for our further analysis in this study. It is noted that BARON requires proper license in order to be deployed.
6 Experiments and Results
6.1 Road Asset Collection
From totally 49,998 vehicle encounters recorded by the University of Michigan Safety Pilot Model Development (SPMD) program, [wenshuo] selected 2,568 naturalistic driving encounters with adequate duration. 10 typical driving scenarios were then successfully extracted. Since there is currently no method that can analytically map traffic encounters to necessary roads and intersections, we extract road assets through manual annotation on OSM-logged road maps. After locating typical driving encounters on OSM via GPS locations, we exported the necessary surrounding road structure in .osm format and imported it into MATLAB using [mathworks:osm] (see Figure 5). Besides, we export the OSM road map of CAV proving ground at Mcity [mcity] and segmented it into 22 road assets which support different typical driving scenarios (see Figure 6). We down-sampled the nodes in OSM raw data to simplify graph representations and assign random integer to value each road asset.
6.2 Solver Performance on Phase 1
We run numerical experiments on BMIBNB, SCIP, BARON to choose the one we will use to proceed. Since BMIBNB is a BNB technique that relies on external solvers, we test it with two solver configurations: a) LP: SCIP, lower bound: SeDuMi, upper bound: fmincon; b) LP: Gurobi, lower bound: Gurobi, upper bound: fmincon. The preliminary testing suite is based on model 4.2.3 with three trivial road asset candidates. We initiate 7 tests with each one challenging the solver to map one non-empty subset of those assets into a given squared space. One single asset is selected in trial 1 to 3; two are involved in trial 4 to 6; and finally all three assets are selected in trial 7. Solving time for each tested solver is recorded and plotted (see Figure 7).
We set the timeout to 200s since the tests are based on trivial road assets with no more than four nodes each. Such timeout turns out to be more than enough because BARON passed every trial in around 1 second. Solving times reaching 200s indicate failed optimization due to timeout. Comparing the performance of two BMIBNB configurations on trial 5 and 7, we notice that the performance of BMIBNB highly depends on the utilized external solvers. SCIP ranks the last with 3 failed trials. We choose BARON to proceed with due to its obvious advantage over alternatives.
6.3 Model Scalability and Complexity Analysis
Using BARON, we explore the scalability of our optimization models with road assets which are extracted from Mcity road map. We gradually increase the number of available road assets until BARON fails to solve within 100s. First, we run feasible mapping searches (phase 1) with 1, 2, and 3 randomly selected road assets. BARON passed all tests with average solving time among 10 trials being 0.1678s, 2.9325s, and 30.9906s respectively for each road asset quantity. Then, we run flexibility optimizations (phase 2) with all three preliminary road assets and most valuable subsets from previous tests. See Figure 8 for example outputs.
Having solved smaller problems successfully, BARON fails in searching for a feasible mapping for more than four assets in phase 1, or in optimizing transition flexibility of more than three assets in phase 2. In both cases, BARON fails to secure any feasible solution during preprocessing. The upper bound consequently stays at default value () after more than 5 hours’ solving. Noticing the significant performance drop in solving, we analyze our model complexity and present the results in this section. We calculate the number of nonlinear constraints and binary variables introduced by each constraint set (see Table 1). Some notations are explained as follow:
: number of nodes in road asset .
/: all internal segments / internal segments in
: cardinality of .
: the number of all transition road candidates.
|Source||Constraints / Binaries||Quantity||Complexity||Type|
|SACS||Dist Cons. 1-3||Quadratic|
|Orient Cons. (5)||Bilinear|
|CACS||Intersection Ineq. 11, 4.2.3||Polynomial|
|Auxiliary Bin. 11||Element-wise|
|Connect.||Intersection Ineq. 4.3.2||Polynomial|
|Auxiliary Bin. 4.3.2-4.3.2||Element-wise|
Analytical result shows that road asset quantity and internal segment quantity are both major factors contributing to model complexity. The usage of excessive binary variables and nonlinear constraints requires a significant amount of branching process and relaxation solving, which has become the major computational burden. Notably, road asset quantity, internal segment quantity, and constraint type are all determined by problem construction and formulation, thereby fixed for a given problem. Therefore, we will work on problem decomposition to address the computational difficulty.
In this study, we proposed the first procedural scenario-based optimization approach for CAV proving ground design, and refer to it as an “Xcity” design problem. The notion of traffic encounters is employed to define road assets as the fundamental construction elements of a CAV testbed. This infrastructure enables the use of scenario-based capability evaluation and direct-connectivity-based flexibility evaluation for CAV proving grounds. Then the Xcity design problem is formulated as general nonlinear optimization problems to fulfill the demanding fidelity requirement of CAV testing and validation tasks. We also presented a two-phase optimization model and demonstrated the design procedures using limited number of road assets. Our numerical results showed the effectiveness of the formulation to guide Xcity constructions strategically. Due to the exponential computational complexity of the resulting formulation, the use of decomposition methods signifies our future direction to extend the tractability of the proposed formulation.
8 Author Contribution Statement
The authors confirm contribution to the paper as follows: study conception and design: Rui Chen, Mansur Arief, Ding Zhao; data collection: Rui Chen; analysis and interpretation of results: Rui Chen, Mansur Arief, Ding Zhao; draft manuscript preparation: Rui Chen, Mansur Arief, Ding Zhao. All authors reviewed the results and approved the final version of the manuscript.