• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    A Distributionally Robust Optimization Method for Passenger Flow Control Strategy and Train Scheduling on an Urban Rail Transit Line

    2022-08-17 07:18:46YahanLuLixingYangKaiYangZiyouGaoHoushengZhouFantingMengJianguoQi
    Engineering 2022年5期

    Yahan Lu, Lixing Yang, Kai Yang, Ziyou Gao, Housheng Zhou, Fanting Meng, Jianguo Qi

    State Key Laboratory of Rail Traffic Control and Safety, Beijing Jiaotong University, Beijing 100044, China

    Keywords:Passenger flow control Train scheduling Distributionally robust optimization Stochastic and dynamic passenger demand Ambiguity set

    ABSTRACT Regular coronavirus disease 2019(COVID-19)epidemic prevention and control have raised new requirements that necessitate operation-strategy innovation in urban rail transit. To alleviate increasingly serious congestion and further reduce the risk of cross-infection, a novel two-stage distributionally robust optimization (DRO) model is explicitly constructed, in which the probability distribution of stochastic scenarios is only partially known in advance. In the proposed model, the mean-conditional value-atrisk (CVaR) criterion is employed to obtain a tradeoff between the expected number of waiting passengers and the risk of congestion on an urban rail transit line.The relationship between the proposed DRO model and the traditional two-stage stochastic programming(SP)model is also depicted.Furthermore,to overcome the obstacle of model solvability resulting from imprecise probability distributions, a discrepancy-based ambiguity set is used to transform the robust counterpart into its computationally tractable form. A hybrid algorithm that combines a local search algorithm with a mixed-integer linear programming(MILP)solver is developed to improve the computational efficiency of large-scale instances.Finally, a series of numerical examples with real-world operation data are executed to validate the proposed approaches.? 2021 THE AUTHORS. Published by Elsevier LTD on behalf of Chinese Academy of Engineering and Higher Education Press Limited Company. This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

    1. Introduction

    The coronavirus disease 2019 (COVID-19) has been greatly affecting people’s daily travel patterns, resulting in significant uncertainty regarding passenger demand in urban public transit.As the main artery of public transit, the urban rail transit system transports a large number of passengers shuttling between different lines and stations every day.In some megacities,the passenger demand in peak hours far exceeds the transport capacity at busy stations (especially at some transfer stations), and a great many passengers must wait on the platforms and are even stranded,resulting in a serious crowdedness issue [1]. According to the World Health Organization, people who leave their home during an epidemic period should refrain from crowding to avoid crossinfection. In other words, passenger congestion in the confined environment of stations necessarily increases the chance of cross-infection and decreases travel comfort. Thus, practical and effective measures—such as limiting the vehicle loading capacity,keeping a safe distance between passengers, and so forth—are needed to ensure safe operation of rail transit with the gradual restoration of social order and the increase in traffic volume brought by the resumption of work and schools.

    At present,two types of measures have been applied in China’s urban rail transit systems, which do indeed mitigate the risk of congestion at stations with large passenger demand: ①The first measure involves passenger flow control. To cope with large passenger demand, the best choice is to impose passenger flow control; that is, to limit the speed of passenger flow entering the platform and make passengers queue in an orderly manner in station halls,so as to avoid passengers gathering on the platforms and thus reduce the risk of cross-infection while improving operation safety.For example,so far,up to 96 stations in the Beijing rail transit system have imposed routine passenger control strategies.②Another practical method is demand-oriented train scheduling;this is a traditional approach used to deal with heavily congested passenger flow in peak hours, which has been studied deeply in the past years[2-4].More precisely,based on detailed information on dynamic passenger demand, this method aims to shorten the train headway to enhance the capacity of the urban rail transit line/system and optimize train schedules so as to minimize the passenger waiting time at stations. However, the phenomenon of huge passenger flow during peak hours has become normative in some megacities, such that even the maximum departure frequency is insufficient to significantly ease the heavy congestion[5]. In addition, as far as most transfer stations are concerned,the number of transfer passengers in practical operations has far exceeded the volume of outside arrival passengers. Xizhimen Station in Beijing can be taken as an illustration. As shown in Fig. 1,the transfer passenger demand at this station is always greater than the outside arrival passenger demand, which is particularly obvious during peak periods. Although some studies [5-9] have paid a great deal of attention to this problem, most still neglect to control transfer passengers simultaneously with outside arrival passengers.To relieve crowding and enhance efficiency as much as possible on an overall line,there is an urgent need to study passenger flow control strategies combined with train scheduling under the influence of transfer passengers—a need that is not addressed in the existing literature.

    In reality,some important information on the urban rail transit system,such as passenger demand,is usually uncertain.To capture this uncertainty, studies [10,11] have used stochastic scenarios with deterministic probability distributions to characterize the randomness of passenger demand, and the scenario-based timedependent passenger demand can be recorded by the automatic fare collection (AFC) system. However, it is difficult to obtain the accurate probability distribution of the involved scenarios in advance, especially in real-world applications [12]. Based on partial known probability information, the modeling methodology corresponds to the distributionally robust optimization (DRO)approach, which has recently been applied to the field of transportation planning and management. Enlightened by this method,the present study aims to find a robust passenger flow control strategy that works in all passenger demand scenarios with partially known probability. For clarity, Fig. 2 shows the passenger control process, in which the time-dependent passenger demand of two scenarios is represented by blue and red curves, respectively, and a robust passenger flow control strategy (represented by the purple line) is generated to control the passenger entering speed and thereby reduce platform congestion. This decisionmaking process is associated with the following potential risks incurred by the two-fold uncertainty: ①Dynamic outside arrival passengers and transfer passengers are both highly uncertain,and are described by stochastic scenario sets; and ②in practice,operators might only know the support of the probability distribution information of the scenarios, while precise distribution information is out of reach. Under these uncertain conditions, risk management [13] has not been sufficiently studied in the existing literature related to the passenger flow control problem,to the best of our knowledge, so it will also be explicitly addressed in our study.

    1.1. Literature review

    Imbalance of supply and demand is the main cause of congestion on an urban rail transit line.Based on this fact,the existing literature has focused on regulating capacity and demand from the supply and demand sides, respectively, in order to alleviate congestion. On the supply side, the literature mainly optimizes the train schedule to meet the dynamic passenger demand with uneven temporal and spatial distribution, so as to avoid a deficiency or insufficiency of transport capacity; on the demand side,a passenger flow strategy,differential pricing,and other measures are mainly used to guide and limit the passenger flow. In addition, the existing studies that focus on dealing with uncertainty in the field of transportation tend to adopt different methodologies, including stochastic programming (SP), robust optimization (RO), and DRO. The following discussion provides a detailed overview of the state of art in the literature.

    1.1.1. Passenger flow control

    Due to the increasing gap between the limited transport capacity and the continuous growth of passenger demand, many urban rail transit lines in China are suffering from overcrowding. The question of how to deal with this problem is thus the focus of operation management departments and many researchers. As a substantial increase in transport capacity is almost impossible in the near future,implementing passenger flow control at congested stations or lines is the best choice to reduce crowdedness. Many researchers have devoted time and effort to this problem, and the related studies can be summarized into the following two categories: The first category is based on the mathematical programming approach and the other one is based on the computer simulation method.

    Fig. 1. Illustration of inbound and transfer passengers at Xizhimen Station.

    Fig. 2. Illustration of the robust passenger flow control process.

    ? Passenger flow control based on the mathematical programming approach.In recent years, an increasing number of researchers have employed the mathematical programming approach to study the passenger flow control problem,which includes three stratification levels: the station level,line level, and network level. Considering the fixed train schedules and the time-dependent passenger demand,existing studies generally aim to minimize the delay time or the waiting time of passengers.For example,Xu et al.[6]proposed a multi-objective optimization model for an urban rail transit station under uncertain demand in order to calculate the number of passengers that should be controlled through the genetic algorithm. In addition, Xu et al. [8] addressed the problem of passenger flow control in an urban rail transit network during peak hours.Their problem was formulated as a mixed-integer programming model to simultaneously adjust the number of inbound and transfer passengers,which could be solved by the genetic algorithm. Shi et al. [9] developed an integer linear programming (ILP) model to generate an optimal passenger flow control strategy for an urban rail transit system. Based on analyses of the alighting and boarding processes, Wang et al. [14] formulated an integer programming model to achieve the optimal state for the entire urban rail transit line,with a focus on minimizing the average passenger delay. Yuan et al. [15] established a mixed-ILP(MILP) model in order to minimize the total passenger waiting time in an urban rail transit network and solved it by CPLEX software. By introducing a timetable-oriented spacetime network representation, Meng et al. [16] proposed an ILP model to produce high-quality passenger flow control strategies for an urban rail transit line.

    ? Passenger flow control based on the computer simulation method.In view of the complex dynamic characteristics of the operation state in an urban rail transit system, computer simulation modeling can reveal the occurrence and evolution of passenger flow congestion, and thus give full play to the advantages of simulation modeling in this problem. Numerous published studies use the simulation method to investigate the passenger flow control at stations. Hoogendoorn and Daamen [17]established microscopic pedestrian dynamics models to assess the design safety in large passenger flow scenarios in Lisbon metro stations.To simulate passenger flow at the Xuanwumen Station (a transfer station) in the Beijing rail transit system during peak hours, Zhang et al. [18] used Vissim software, which achieved good performance. Seriani and Fernandez [19] studied the effects of pedestrian traffic management on passengers’ boarding and alighting time,and used experimental means to obtain the standards for pedestrian traffic management at the platform and doors of metro vehicles.Fei and Liu[20]used Anylogic software to formulate a simulation model of passenger flow organization and accurately identified the congestion points according to the simulation results. Jiang et al. [21] proposed a coordination model of passenger flow control and train jumping based on the Q-learning approach,which could be efficiently solved by simulation methods.

    1.1.2. Train scheduling in urban rail transit

    Over the past several decades, a large and growing body of literature has proposed different mathematical models and solution algorithms for train scheduling problem in urban rail transit. For example, to minimize passenger waiting time, Niu and Zhou [4]developed a nonlinear programming model to optimize the train schedule for the whole line, which could be optimized by the genetic algorithm. They then took oversaturated situations into consideration and formulated a mixed-integer nonlinear programming model[22].Li et al.[23]studied the fairness issue in the train scheduling process. To optimize train schedules by considering min-max fairness and α fairness, they established a mixedinteger programming model that could be solved by a simulated annealing-based adaptive large neighborhood search algorithm.Yin et al.[24]proposed an MILP model to fix the coordinated train scheduling problem with unlimited train boarding capacity. Yang et al. [25] developed a 0-1 ILP model for the last train timetable in an urban rail transit network, which could be effectively solved by a heuristic algorithm based on Lagrangean relaxation. Liu et al.[26] formulated a mathematical optimization model to determine a robust train schedule for maximizing robustness and minimizing total energy consumption, so as to avoid delay propagation as much as possible. Huang et al. [27] considered the rescheduling problem in urban rail transit systems. Two mixed-integer nonlinear programming models with different recovery strategies were formulated to reschedule trains during track disruption, which could be linearized by a hybrid approach combining big-M and time-indexed formulations. Tian and Niu [28] addressed the demand-oriented train scheduling problem under overtaking operations. Their problem was formulated into an ILP model, and a corresponding novel surrogate-dual-variable column generation algorithm was proposed to generate approximate optimal solutions.Mo et al. [29]proposed a mixed-integer nonlinear programming model to generate an optimized train timetable with suitable train speed profiles, as well as an efficient rolling stock plan, aiming to minimize both passenger waiting and travel times. Qi et al.[30]formulated an ILP model for the integrated stop planning and timetabling problem, with the aim of minimizing the total travel time of trains.

    Overall, the above studies address the passenger flow control problem and the train scheduling problem separately, thus failing to achieve an optimal system effect or reduce the pressure on urban rail transit operation to the greatest extent. It should be noted that integrated optimization of the above two aspects is much more challenging than dealing with them separately because of the complex coupling relationship between trains and passenger flow.To date,only three papers can be found in the literature that use a mathematical programming approach to carry out an integrated study of these two aspects: Xu [31] formulated a quadratic programming model to generate a multi-station collaborative passenger flow control strategy and a corresponding train schedule;in this model, the decision variables include the number of inbound passengers and the dwell time of trains at certain key stations.Shi et al.[5]established an ILP model to obtain an optimal passenger flow control strategy for the whole line and a corresponding train schedule to improve the service quality,in which the passenger flow control strategy was implemented on each timestamp.Gong et al. [32] implemented the passenger flow control based on the train status. However, none of these studies consider the influence of transfer passengers.

    1.1.3. Approaches for dealing with uncertainty

    It is well known that transportation activities usually faces a lot of dynamics and uncertainties due to system complexity[1,11,33,34]. With the development of optimization theory and the improvement of computing technologies, an increasing amount of attention has been paid to how to deal with uncertain parameters. Thus far, three major kinds of approaches—namely,SP, RO, and DRO—have been developed to handle uncertainty.More specifically,the SP method assumes that the probability distribution of uncertain parameters is known in advance. For example, Gong et al. [35] considered the stochastic environment in an urban rail transit system and developed an integer nonlinear programming model to simultaneously optimize the number of train services, the headway settings, and the speed profile selection decisions. Errico et al. [36] constructed a two-stage SP model for the vehicle routing problem with hard time windows,which could be solved by exact branch-cut-and-price algorithms. It is widely recognized that the SP method is based on a given probability distribution.However,the real value of probability distribution is difficult to obtain in practice and usually needs to be estimated through historical data. Due to the existence of estimated errors,major deviations may occur in practical applications, which is a great weakness of this method. The RO approach considers that uncertain parameters are all based on uncertainty sets without any distribution assumption. The aim of the RO approach is to overcome the uncertainty and find the best possible solutions that satisfy all constraints. For example, to handle the uncertainty of passenger demand, Yang et al. [10] proposed an MILP model with the max-min reliability criterion, aiming to increase the number of last-train successful transfer passengers and reduce the total running time. Qi et al. [37] studied the train timetabling and stop-planning problem with uncertain passenger demand. An ILP model was formulated and solved by CPLEX software. Based on the technique of Light Robustness, Cacchiani et al. [38] proposed three MILP models to derive robust train-stop plans and timetables, in order to reduce passengers’ inconvenience. Here, we note that a serious drawback of the RO approach is that it is too rigid and conservative; however, it does not require the distribution information to be known in advance, unlike the SP method.

    The rapid development of the DRO approach has overcome the aforementioned shortcomings of the SP and RO methods.The DRO approach combines the statistical learning technique with optimization theory;it can obtain a good enough solution by assuming that the parameters obey some possible distributions.The theories underpinning the DRO approach have been well studied over the last decades [39,40]. For example, Delage and Ye [39] produced one of the first studies on this approach,in which the DRO models were formulated under the first-order and second-order moment inequalities, respectively. Esfahani and Kuhn [40] proved that the dual form of a Wasserstein-metric DRO model could be decomposed into finitely convex optimization problems with tractable computational complexity. Recently, this approach has also attracted much attention from certain engineering fields [41-47].For example, focusing on the planning and scheduling problem under uncertain demand, Shang and You [41] formulated a twostage DRO model to produce less conservative solutions. For finite and infinite horizon optimal control problems, Van et al. [42] proposed distributionally robust constrained formulations to provide frameworks that could generate robust control policies. Considering the integration of urban and rural public transport systems,Shang et al. [43] developed a DRO model with ambiguous chance constraints regarding the travel time requirement, with the aim of minimizing the total construction cost.For a disaster relief management problem under demand uncertainty,Wang et al.[44]proposed a DRO model to simultaneously optimize the integrated facility location, inventory pre-positioning, and delivery decisions.Yang et al. [45] presented a distributionally robust chanceconstrained program to model the last-train coordination planning problem in an urban rail transit system,in which uncertain passenger demand was captured by an ambiguous set. Wang et al. [46]focused on the hub location problem with multiple commodities under uncertain demand and cost. This problem was modeled as a DRO model, which was further reformulated as a moderatesized MILP for ease of calculation. Zhang et al. [47] established a DRO model with the Wasserstein distance-based ambiguity set for the vehicle routing problem with time windows, aiming to minimize the risk of late customer arrivals.

    Thus far, situations with uncertain passenger demand have received less attention in the area of passenger flow control. As far as we know, there is only one relevant reference in this field that addresses uncertain demand: Meng et al. [11] proposed an ILP model on an urban rail transit line by applying the SP method;in this model, uncertain passenger demand was described using a stochastic scenario set with completely pre-given distribution information. Moreover, there is no relevant theoretical achievement employing the DRO approach in the field of passenger flow control combined with train scheduling.

    In summary,the passenger flow control problem has been a hot research topic in recent years,and some mathematical models and simulation methods to address this problem have been well studied. Nevertheless, most of these approaches still separately consider or partially integrate the three aspects of passenger flow control, train scheduling, and SP method. Moreover, no existing study simultaneously deals with the two-fold uncertainty with respect to the passenger demand and scenario probability, even though it is a significant and challenging issue for real-world applications, as it is difficult to predict passenger demand accurately.With this concern,this paper is particularly interested in developing a DRO method for the passenger flow control strategy and train scheduling on an urban rail transit line,as well as a powerful algorithm with fast computing speed.

    1.2. The contributions of this study

    To fill the abovementioned gaps, this paper aims to make the following contributions to the study of passenger flow control and train scheduling problem on an urban rail transit line:

    (1)A DRO model,which focuses not only on both outside arrival passengers and transfer passengers, but also on train schedules, is rigorously formulated. Compared with the existing literature, this model has two significant improvements: ①The proposed model considers the inherent uncertainty of passenger demand. Since passenger demand in the real world is usually uncertain and dynamic, a series of stochastic scenarios are developed to depict this uncertain parameter in this paper. ②The precise probability distribution information of each stochastic scenario does not need to be known completely in advance—that is,only the partial distribution information is needed, which is the main reason why this approach is employed in this paper.Since very few existing studies consider the inherent uncertainty of passenger demand and the characteristics of scenario probability cannot be predicted in advance, this is a novel idea.

    (2) Our distributionally robust formulation is a semidefinite programming model, which is a notoriously difficult problem to solve.To convert this model into a computationally tractable form,a discrepancy-based ambiguity set is employed. Furthermore, the equivalent MILP model can be obtained through the strict theoretical proofs. Based on an analysis of the complexity of our formulation, we propose an efficient hybrid algorithm by combining a local search algorithm(LS)with a suitable MILP solver to generate high-quality solutions in an acceptable computation time, especially for large-scale instances.

    (3) Finally, to demonstrate the effectiveness and practicability of our proposed approaches, a series of numerical experiments are implemented based on operating data from Line 15 in the Beijing rail transit system. We compare the experimental results of our DRO model with those of the traditional SP model and the classic RO model.These experiments demonstrate the performance and advantages of the proposed optimization model and solution approaches.

    The remainder of this paper is organized as follows.In Section 2,we give a detailed description of the passenger flow control and train scheduling problem on an urban rail transit line.In Section 3,we formulate the problem into a two-stage SP model and a DRO method, respectively. Then, Section 4 proposes a hybrid algorithm based on the LS and a suitable MILP solver.In Section 5,we implement a case study,that is,a large-scale case based on the operating data of Line 15 in the Beijing rail transit system,to show the superiority of our proposed approaches. Finally, Section 6 gives some conclusions and future research.

    2. Problem description

    Consider an oversaturated urban rail transit line,which consists of a series of common and transfer stations(Sorand Strare used to denote the sets of common and transfer stations,respectively),and a set of physical segments between stations (the kth segment is defined as the connection from stations k to (k + 1)). Denote the set of stations on this line as S, then S=Sor∪Str= 1,2,..., |S |{ }. In operations, the transfer between different lines is a common passenger behavior at each transfer station. In order to ease congestion, we aim to collaboratively optimize the train schedule and passenger flow control strategy on this oversaturated line by taking the transfer passenger flow into account.For descriptive clarity,an illustrative urban rail transit system is depicted in Fig.3,where Line 1 is the focus line, and Lines 2 and 3 are the adjacent lines.Typically, Stations 2,3,5 ∈Sor, and Stations 1,4 ∈Str. Passengers on Lines 2 and 3 are allowed to transfer to Line 1 at Stations 1 and 4,respectively.At Stations 2,3,and 5,the inbound passengers are mainly from the outside of the transit system.At Stations 1 and 4, both outside arrival passengers and transfer passengers from Lines 2 and 3 should be taken into consideration, in the process of optimizing the passenger flow control strategies.

    For an urban rail transit line with uncertain passenger demand,congestion is directly related to the following four factors:outside arrival passengers, transfer passengers from adjacent lines, the schedules of involved trains,and the train capacity.Due to the limited improvement of train capacity in real-world operations, the most effective method to alleviate congestion is to consider the effects of the other three factors from the system perspective,and to carry out collaborative optimization of the train schedules and passenger flow control strategies. In the process of passenger flow control, outside arrival passengers are first required to queue in station halls,and they are then released to the platform to board train i under the robust passenger flow control strategy at station k in scenario ω (denoted by(ω )). Nevertheless, the transfer passengers from adjacent lines (i.e.,(t ,ω)) are allowed to go to the platform freely and wait for the next train arrivals, as shown in Fig. 4. Once the number of outside arrival and transfer passengers exceed the capacity of the trains, outside arrival passengers who come later than this time must be detained in the station hall; they then gradually enter the platform according to the robust passenger flow control strategy. In this paper, the control strategy is only imposed on outside arrival passengers,while transfer passengers are released to the platform and can take trains freely. In addition, all the passengers on the platform are required to take the train that arrives earliest, without waiting for the next train, in order to guarantee the operation safety.

    In summary, the focus of this paper is to develop a multiscenario mathematical formulation for the robust passenger flow control problem combined with train scheduling on an urban rail transit line,with the goal of minimizing the sum of operating time,the expected number of waiting passengers, and the conditional value-at-risk (CVaR) of congestion. Our formulation also explicitlycaptures the time-dependent and stochastic passenger demand—that is, the outside arrival passengers and the transfer passengers—and the train capacity, aiming to achieve the optimization of system performance. The proposed mathematical formulation is based on the following three assumptions.

    Fig. 3. An illustrative urban rail transit system.

    Assumption 1:The outside arrival and transfer passenger demand is known in advance by processing the historical AFC data.A similar assumption has been widely adopted in previous studies[1,8,11].

    Assumption 2:The section running time and station dwell time of trains are pre-given by rail managers according to the actual operating data. That is, these parameters are not affected by the numbers of boarding or alighting passengers.

    Assumption 3:We only consider the direction from station 1 to station |S |as the operation direction,while the other direction can be operated according to the same optimization method. In addition, the inbound and transfer passengers, who are released to platforms under the robust passenger flow control strategy, can all be taken away by the first arrival train.This assumption is common in passenger flow control and train timetabling studies[4,32].

    3. Mathematical formulations

    In this section,we first introduce a two-stage SP model for passenger flow control and train scheduling problem with completely known probability information in Section 3.1.We then formulate a novel DRO model with uncertain probability in Section 3.2.

    3.1. Two-stage SP model with completely known probability information

    To characterize the above problem in a stochastic version, the modeling process will be explicitly discussed in the following sections. In Section 3.1.1, we first list all the related notations and decision variables in the formulation.The events that describe the dynamic processes of passenger loading and train operations are introduced in Section 3.1.2. In Section 3.1.3, we present a two-stage SP model based on the mean-CVaR criterion.

    3.1.1. Notations and decision variables

    In this paper,we use boldface letters to represent vectors and V is employed to represent the set of random variables. Uncertainty is modeled via a set of stochastic scenarios, the probability distribution of which is denoted as p[·]. We use Ep[x ] to represent the expectation of x with probability distribution p[x]. Based on these definitions,we next develop a two-stage SP model for the problem of interest, in which the first stage determines the train schedules and the second stage is associated with the movement of passengers. To formulate the model rigorously, five sets of decision variables are defined as follows.

    Decision variables:

    Sets:

    Parameters and notations in the first-stage model:

    Parameters and notations in the second-stage model:

    3.1.2. Passenger loading and train operations

    For the sake of clarity, this section introduces a non-increasing binary variable di,k(t) to denote the departure characteristics of train i at station k and timestamp t.Specifically,di,k(t)=0 indicates that train i has departed from station k at timestamp t,and di,k(t)=1 indicates otherwise.As shown in Fig.5,train i departs from station k at timestamp 2. It is clear that a series of ‘‘1” for di+1,k(t) - di,k(t)indicates the departure headway between trains i and(i+1)at station k,during which passengers can board train(i+1)according to the robust passenger flow control strategy. By using this binary variable, it is possible to formulate the problem of interest as an optimization model with linear forms, as described in the following discussion.

    3.1.3. A two-stage SP model

    A two-stage SP model based on the mean-CVaR criterion is formulated in this section, in which the distribution of the uncertain passenger demand is assumed to be completely known in advance and is captured by multiple scenarios with deterministic probability distribution.In the first stage,we propose a formulation that is only relevant to train scheduling, shown in the equations below:

    Fig. 5. Illustration of scheduling variables at station k.

    where Q (d ,ξ (ω )) is the random optimal value of the second-stage model, d denotes the first-stage decision vector of train schedules,ξ (ω ) represents the random input data of passenger flows, E is the expected number of waiting passengers,λ is the tradeoff coefficient that reflects the decision-makers’ risk preference, ζ1and ζ2are two cost coefficients, and CVaR denotes the conditional value-at-risk(VaR) criterion, which is a popular risk measure to model the decision-maker’s risk preference. Compared with the VaR criterion,CVaR highlights the average level of excess losses that occur beyond the VaR cutoff point in the distribution.Therefore,CVaR is superior to VaR in reducing uncertainty risk,since more risk-averse decisions can be obtained by adopting the CVaR criterion.This could be particularly helpful in preventing serious situations in which a large number of passengers congregate at an urban rail transit station.

    The objective function in Eq. (1) aims to minimize the sum of the operating time (i.e., T0), the expected number of waiting passengers, and the CVaR for making decisions. In particular, if λ = 0,the obtained solution will be the result of not taking the CVaR for making decisions into consideration. It is clear that this objective function will minimize the operating time if ζ2is set to zero.In this case,this problem turns out to be a supply-side-oriented optimization. However, when considering multi-scenario coupling,only minimizing the operating time will inevitably increase the risk of service unfulfillment. Thus, the CVaR criterion is employed to handle the decision risk caused by uncertainty,and the number of waiting passengers is the most common decision criterion in the passenger flow control problem.

    Given the section running timeand the station dwell time,Eqs.(2)and(3)formulate the arrival and departure times for each train.In particular,is an input parameter that is fixed in advance.Eq. (4) defines the departure headway between two consecutive trains i and (i + 1) arriving at station k. Eq. (5) guarantees the minimum and maximum safe departure headway. By imposing Eq.(6),we can obtain the actual time with respect to the departure headway.Eq.(7)ensures that the 0-1 decision variable di,k(t)is nonincreasing for train i at station k.Eq.(8)is employed to calculate the total operating time for the focus line by summing up the headway,station dwell time,and section running time of each train at all stations.Finally,the domain of the decision variable,di,k(t),is defined by Eq.(9).

    In the second stage, the corresponding robust passenger flow control strategy is optimized, which is formulated as follows:

    In this stage,we need to formulate the relationship between the train schedule and the robust passenger flow control strategy, by establishing their coupling constraints. The objective function Eq. (10) aims to minimize the total number of waiting passengers for different scenarios.As we consider the direction from station 1 to station |S| as the operation direction, and thus station |S| is the terminal station with no passenger demand,the number of waiting people at this station is zero. By imposing Eq. (11), it is ensured that all passenger demand should be satisfied over the considered time horizon T.The number of passengers waiting for train i in the station hall can be calculated by Eq. (12). Since all outside arrival passengers need to queue in the station hall and wait for boarding permission, the number of waiting passengers is the difference between the cumulative number of outside arrival passengers and the cumulative number of passengers being served. Eq. (13)tracks the number of transfer passengers boarding train i at station k.As the transfer passengers can directly reach the platform,these passengers must have a high priority to get on trains. By using Eqs. (12) and (13), passenger demand and train schedules can be tightly connected in a series of unified formulations. Eqs. (14)and (15) require the number of waiting passengers to be nonnegative because the number of passengers who are allowed to get on any train i cannot be greater than the number of outside arrival and transfer passengers, respectively.

    Next,we analyze the dynamic passenger loading process,which consists of two types of activities: boarding and alighting. Since multiple stochastic scenarios are considered in this study,a robust passenger flow control strategy,(ω ), which is suitable for all involved scenarios,would be generated.However,due to the diversity of passenger demand, an actual control strategy,(ω ), also exists in each respective scenario. Eq. (16) is used to restrict the number of outside arrival passengers actually boarding train i,which should not exceed the number of waiting passengers.Through Eq.(17),the total number of passengers actually boarding train i at station k can be calculated. Eq. (18) is imposed to ensure that the actual number of boarding passengers is less than or equal to that value under the robust passenger flow control strategy.Eq. (19) ensures that the number of outside arrival passengers allowed to board train i under the robust passenger flow control strategy is equal for each scenario, which greatly improves the chance of practical application in the real world. Eq. (20) makes it possible to obtain the total number of passengers that are allowed to board train i at station k under the robust passenger flow control strategy.Specifically,at common stations,the boarding passengers under the robust passenger flow control strategy only include outside arrival passengers,while the transfer passenger demand is also required to be considered at transfer stations. Eq. (21) makes sure that the number of alighting passengers under the robust passenger flow control strategy from train i at station k is equal to the number of passengers heading to this station, who are in train i in front of this station. Eq. (22) considers the dynamic variation of in-vehicle passengers under the robust passenger flow control strategy when train i leaves station k.More specifically,when train i arrives at its start terminal (station 1) with empty vehicles, the number of invehicle passengers will be equal to the number of boarding passengers under the robust passenger flow control strategy.This value is equal to 0 when train i leaves the terminal station because all passengers are required to get off the train. In comparison, at each intermediate station,the number of in-vehicle passengers is closely related to the total number of in-vehicle passengers at the previous station, the number of alighting passengers at the current station,and the number of boarding passengers at the current station.Moreover, the number of in-vehicle passengers cannot exceed the capacity of the trains,as expressed in Eq.(23).

    If the discrete probability distribution of different scenarios can be specified, we can formulate the expected number of waiting passengers as a computationally tractable form, which can be expressed as follows:

    In practice,it is rather difficult to precisely obtain the probability distribution of stochastic scenarios.In such a circumstance,the probability distribution is assumed to be only partially known; in this study, it belongs to an ambiguity set, which will be analyzed in the following sections.

    3.2. DRO model with uncertain probability

    Based on the notations in the previous sections, we formally propose the following DRO model in Section 3.2.1; next, in Section 3.2.2, we specifically design a discrepancy-based ambiguity set for our developed model. In Section 3.2.3, we further reformulate the DRO model as an MILP model.Finally,we analyze the complexity in terms of the number of decision variables and key constraints in Section 3.2.4.

    3.2.1. DRO model

    Given the uncertain probability of each stochastic scenario, a DRO model based on the worst-case mean-CVaR criterion is formulated in this section, as shown below:

    Below, we discuss the equivalent formulations for the DRO model. First, the objective function in model (Eq. (28)) can be transformed into

    Note that the DRO model (Eq. (28)) is a semidefinite programming, which is computationally intractable for general ambiguity sets. Hence, to solve it directly using a commercial solver (e.g.,CPLEX) within an acceptable computation time, we can transform the DRO model (Eq. (28)) into its computationally tractable form.In the next section, a special ambiguity set is designed to characterize the distribution uncertainty, and the corresponding computationally tractable forms of the DRO model(Eq.(28))are deduced.

    3.2.2. Ambiguity set

    In this study, uncertainty is modeled via a series of stochastic scenarios, and the corresponding probability distributions are characterized by a discrepancy-based ambiguity set [49], that is,P, which is defined as follows:

    To be specific,p0is a vector representing the nominal distribution of the discrete probability, ? denotes the fluctuation vector,and Ψ is a real value between 0 and 1. Note that the condition eT?=0 is formulated to ensure that the vector p meets the requirements eT(p0+?) = 1. That is, the discrepancy-based ambiguity set describes the deviation from a reference(often empirical)distribution. The reasons why we use this ambiguity set are threefold:①This set is simple and practical. The former part, p0, can be obtained from the data of different dates in practical operations,which is easy to realize in engineering. ②This set is convenient for transformation into a computationally tractable form by applying the related theoretical analysis. ③This set often possesses asymptotic properties [50].

    3.2.3. Deterministic equivalent formulation

    Under the designed ambiguity set,the deterministic equivalent formulation can be obtained by invoking the theoretical results of Ref. [49], which are shown as follows.

    Theorem 1.Under the discrepancy-based ambiguity set P, the equivalent form of the DRO model (Eq. (28)) can be expressed as follows:

    where μ,η,γ,μ′,η′,γ′∈R×R|W|×R|W|×R×R|W|×R|W|are auxiliary variables or vectors, and Ψ = eΨ.

    Proof. Under the discrepancy-based ambiguity set P, Eq. (32)can be transformed into the following equivalent form:

    With the strong duality theory of linear programming, the dual form of Eq. (32) can be formulated as the following MILP model:

    Likewise, the dual programming of Eq. (33) is developed as the following linear programming model:

    By regrouping the terms,model(Eq.(39))can be represented as model (Eq. (35)). The proof of Theorem 1 is thus complete.

    Remark 2:In the case of Ψ = 0, the DRO model (Eq. (35)) is equivalent to the two-stage SP model (Eq. (27)).

    3.2.4. Complexity analyses

    The decision variables in the proposed DRO model(Eq.(35))can be categorized into three categories: The first category comprises the variables associated with the train schedule—that is,,,hi,k,di,k( t)—and the second is related to the robust passenger flow control strategy—that is,(ω ) and( ω). The last category comprises the auxiliary variables, such as μ, η(ω), γ(ω),μ′, η′(ω), and γ′(ω). Obviously, the scales of the abovementioned decision variables are fully dependent on the number of stations,trains, timestamps, and scenarios. We are particularly interested in analyzing the complexity of the proposed model (Eq. (35)) in detail,as shown in Table 1.To illustrate this problem more clearly,an instance with 40 in-service trains,20 stations,180 timestamps,and three scenarios is provided in order to identify the total number of decision variables and constraints.In this case,there are over 150 000 decision variables and 319 980 constraints, which indicates that the proposed model (Eq. (35)) is a fairly complex problem.

    4. Solution approaches

    Based on the analyses above,if the scale of the problem is small,the proposed model(Eq.(35))can be solved by means of MILP solvers(e.g.,CPLEX)in a short time,due to the linearity of the model.However, it is difficult to solve real-world large-scale instances.Thus,we need to design heuristic algorithms to search for approximate optimal solutions within an acceptable computation time,in order to meet the needs of real-world applications. Therefore, an effective hybrid algorithm combining an LS with an MILP solverwill be developed in this section to find high-quality solutions.The framework and the detailed techniques of our hybrid algorithm are introduced in the following discussion.

    Table 1 Numbers of decision variables and constraints in the DRO model.

    4.1. Algorithm framework

    An algorithm framework is developed to solve the proposed DRO model, following the procedure shown in Fig. 6. In brief, the original problem (Eq. (35)) is decomposed into two subproblems.The first subproblem (denoted by MP) determines the decision variables related to the train schedule for the focus line, while the second subproblem (denoted by SCP) optimizes the control variables for uncertain passenger demand over all the involved scenarios. An LS is designed to search for feasible train schedules.Then, each feasible schedule is taken as the input for subproblem SCP to produce the corresponding robust passenger flow control strategy, which is used to evaluate the generated schedule in subproblem MP. We observe that subproblem SCP is an MILP model that can theoretically be solved to optimality. However, similar to the problem in Ref. [11], the scenarios in subproblem SCP are tightly coupled across different scenarios, making it difficult for the problem to be solved directly by MILP solvers. Therefore, we decompose it here into a series of single-scenario-based subproblems SSP in order to speed up the computation, in which the train schedule is the input for subproblem SCP and SSP. As each subproblem SSP is an MILP problem and does not involve the coupling relationship among different scenarios, it can be efficiently solved using a suitable MILP solver.

    Next, the minimum value of the optimal solutions to subproblems SSP, denoted by, is treated as the input to subproblem SCP. That is, we first let

    The solution space of the decision variableω( ), which is the key variable that couples the various scenarios together, is greatly reduced by Eq.(41);thus,Eq.(41)acts as a so-called‘‘valid inequality”to the problem of interest.The proposed valid inequality is then taken as the input for subproblem SCP.

    4.2. Local search algorithm

    Fig. 6. The procedure for solving the DRO model.

    To solve subproblem MP,an initial solution should first be generated. Next, the approximate optimal solution is searched in a greedy way. During the searching process, only the best solution in the neighborhood space can be accepted as a seed solution.Finally,the algorithm is terminated when the current best solution cannot be updated within a given maximum number of iterations,or when the number of search iterations reaches a predetermined threshold. Then, the best solution that has been encountered is outputted as a near-optimal solution. For the sake of description,we introduce some important notations in the algorithm: n denotes the index of iteration, nmaxdenotes the threshold of the number of iterations, M denotes the maximum number of iterations for which the encountered best solution is not updated, Xndenotes the best solution at iteration n,X0denotes the initial solution, m denotes the number of candidate solutions at each iteration, N(Xn-1) denotes the neighborhood at iteration n, X* denotes the current best solution, f* denotes the current best objective value, and fndenotes the best objective value at iteration n.

    4.2.1. Initial solution generation

    In this study, the initial solution is generated randomly as the seed solution for the first iteration.That is,we first generate a total of(|I |ˉ1)headway variables randomly within a given rangeand thus an initial solution X0can be produced, whose feasibility can be tested by Algorithm 1, shown below. If this solution is feasible, it will be set as the initial solution; otherwise, we need to regenerate the headway vector until a feasible solution is found.In the searching process, the initial solution will be treated as the input data for subproblem SCP, and the corresponding objective value is used as the evaluation of this train schedule.

    4.2.2. Neighborhood generation

    Algorithm 1. The checking procedure for solution feasibility.

    It is worth mentioning that we need to test the feasibility of each newly generated neighbor solution by means of Algorithm 1.If the solution is feasible,we then insert it into the neighborhood of the current seed solution; otherwise, a new solution must be regenerated until its feasibility is guaranteed.

    4.2.3. Algorithm procedure and termination criteria

    The procedure of the LS + MILP solver algorithm is reported below as Algorithm 2, which is designed based on the aforementioned specific techniques.Two criteria are employed to terminate the searching process:①If the current iteration n is greater than a pre-given threshold nmax,the search process should be terminated.②If the objective value of the currently encountered best solution is not updated within M iterations, the LS is terminated, and the current best solution can be treated as the near-optimal solution.We use y to denote the number of iterations that the current best solution is not updated.

    5. Numerical experiments

    In this section,a series of numerical experiments based on realworld operation data from the Beijing rail transit system are conducted to demonstrate the effectiveness and performance of ourproposed models and algorithm. More specifically, the stochastic and dynamic passenger demand is obtained by processing historical AFC data,and the running time of each section is deduced from the practical timetable.To solve this large-scale problem within an acceptable computation time, we use the proposed algorithm framework to generate near-optimal solutions for the models.The algorithm is coded by MATLAB 2014a with CPLEX 12.6 on a Windows 10 personal computer with Intel Core i5-10500 central processing unit (CPU) and 16G random-access memory (RAM).

    Algorithm 2. The procedure of the LS + MILP solver algorithm.

    5.1. Data description

    As shown in Fig. 7, this set of experiments consider Line 15 of the Beijing rail transit system(i.e.,the blue line),which has a total length of 41.4 km and 20 stations. This line connects a residential zone in a suburban town with the downtown working area, and numerous commuters need to ride this line on weekdays,resulting in over-saturation in the morning peak hours in the downtown direction. We only consider the direction from station Fengbo(FB) in the suburbs to station Qinghuadongluxikou (QH) downtown as the experiment environment, for which the station dwell time and the section running time are taken from real operation data, as shown in Table S1 in Appendix A. We select the morning peak hours from 7:00 to 10:00 as the considered time horizon.To characterize the dynamics of the passenger demand, the discretized time granularity is set as tunit= 60 s in this series of experiments, and a total of 180 timestamps are involved (i.e.,T = {1,2,...,180}),so as to balance the real-world operation conditions and the computation efficiency. To depict the uncertainty of passenger demand, three typical stochastic scenarios with unknown probabilities are taken into account, for which the nominal probability distribution is predetermined to be=(0.20,0.30,0.50).Todeterminetheambiguity set, thevalue of the adjustableparameterΨis setas0.02,0.06,and 0.10,respectively; the risk-aversion parameter α is chosen from the set{0.95,0.50,0.05}; and the tradeoff coefficient λ is chosen from the set {0.10,0.50,0.90}. For illustrative clarity, the dynamic outside arrival passenger demand in scenario 1 is displayed in Fig. 8.In addition,it should be noted that there are four transfer stations on this focus line:station Wangjing(WJ)connects Lines 15 and 14,station Wangjing West (WJX) connects Lines 15 and 13, station Datunlu East(DTLD)connects Lines 15 and 5, and station Olympic Green (ALPK) connects Lines 15 and 8. These four stations attract many transfer passengers in addition to the outside arrival demand. In the data preparation stage, the passenger ratios for the outside arrival and transfer passengers at different destinations are given in Tables S2 and S3 in Appendix A, respectively. In the algorithm, the threshold of the number of iterations nmaxis set as 100.The other relevant parameters are listed in Table S4 in Appendix A,where C is set as the capacity using the 120%full-loading rate of each in-service train.

    Using the above parameter settings, a series of numerical experiments are conducted to demonstrate the performance of our proposed approaches. The first set of experiments are carried out to evaluate the validity of the proposed DRO model in Section 5.2. The experiments in Section 5.3 aim to compare the performance of the SP model(Eq.(27))with that of the DRO model(Eq. (35)). In addition, the experiments in Section 5.4 are used to evaluate the computation results of the classic RO model(Eq. (S2) in Appendix A) and the DRO model (Eq. (35)). Finally,Section 5.5 reports the impact of some key coefficients in the final set of experiments.

    5.2. Computation results of the DRO model

    To test the proposed DRO model (Eq. (35)), we consider the computation results by using different combinations of parameters α,Ψ,and λ.It should be noted that the cost coefficients ζ1and ζ2in the objective function are used to make tradeoffs between operations and services. Here, to make the comparison of optimization results more obvious and to consider the difference in the order of magnitude for two objectives, we set ζ1as 100 000 and ζ2as 1 in this set of experiments, and the detailed sensitivity analyses related to these two coefficients are shown in Section 5.5.

    The computation results are displayed in Table 2. For clarity,Fig. 9 presents the optimal objectives obtained by the DRO model with different combinations of parameters. Typically, with the increase of parameter Ψ, the optimal objective of the model(Eq. (35)) gradually worsens. Furthermore, the tradeoff coefficient λ and the optimal objectives have the same monotonicity, and the optimal objective increases when the risk-aversion parameter α increases. These phenomena are consistent with the theoretical results. More specifically, the parameter α represents the level of risk aversion;that is,α=1 represents extreme risk aversion.Therefore, a larger parameter α leads to a more conservative decision plan,which corresponds to a larger objective value.In the ambiguity set,Ψ represents the disturbance range of the uncertain probability. As the parameter Ψ increases, the distribution of probability information to the decision-maker becomes increasingly coarse, resulting in much worse outcomes for the decision.

    Fig.7. Illustration of Line 15 of the Beijing rail transit system.QH:Qinghuadongluxikou;LDK:Liudaokou;BST:Beishatan;ALPK:Olympic Green;ALL:Anlilu;DTLD:Datunlu East; GZH: Guanzhuang; WJX: Wangjing West; WJ: Wangjing; WJD: Wangjing East; CGZ: Cuigezhuang; MQY: Maquanying; SH: Sunhe; GZ: Guozhan; HLK: Hualikan;HSY: Houshayu; NFX: Nanfaxin; SM: Shimen; SY: Shunyi; FB: Fengbo.

    Fig. 8. Illustration of the outside arrival passenger demand of Line 15.

    For clarity,Fig.10 displays an optimized train schedule and the number of in-vehicle passengers in each train over each section,where the red lines represent the full-load sections. Clearly, most departure headway during the 7:00-7:42 time period is relatively small,and the full-load situations are all associated with the trains that depart from the start terminal station during this time period,which demonstrates that the train schedule and passenger demand are closely coupled. Next, we are particularly interested in investigating how the objective value converges to the optimal value in the searching process of our proposed algorithm. Since the convergence tendency of the objective values is similar for different instances,we only present the results of the instance with α=0.95,λ=0.10,and Ψ=0.02 for illustrative convenience here,as shown in Fig. 11.It follows from this figure that the best objective value can be iteratively reduced by the proposed LS + MILP solver algorithm,where the solution quality is gradually improved with a fast convergent speed. More specifically, the best objective value decreases drastically in the first several iterations, but is not updated any longer in the last few iterations. Thus, we conclude that the LS + MILP solver algorithm exhibits a good performance in solving the developed DRO model(Eq.(35))in terms of solution quality, while requiring a reasonable computation time.Furthermore, to clarify the impact of the tradeoff coefficient, λ, in the following context, we will discuss the variations of the expected number of waiting passengers and CVaR when parameter λ changes. We summarize this section with the following observation.

    Observation 1.(The necessity of considering different combinations of α, λ, and Ψ.) It follows from the computation results that the combination of parameters α, λ, and Ψ has significant effects on the DRO model.If the decision-makers prefer a lower risk decision,it is better to set a bigger α and λ,and to know the probability distribution as precisely as possible (in which case, Ψ should be small).That is,the level of risk aversion and the disturbance range of uncertain probability are significant factors in the process of optimizing the train schedules and passenger flow control strategies.

    Next,we further consider the variation of the expected number of waiting passengers and CVaR with respect to different values of parameter λ.For comparison,we set parameter α as 0.50 and set Ψ as 0.02 and 0.06,respectively.Fig.12 depicts the variation trend of the expected number of waiting passengers and CVaR with different values of λ. Clearly, for two instances in this figure, the CVaR yields a downward linear increasing tendency in the optimal objective if we increase the coefficient λ, while the expected value shows the opposite trend when α and Ψ are fixed. In addition,when λ is set as 0.5, the risk preference of the decision-makers is neutral. Thus, the expected value and CVaR are very close in thiscase,as shown in Fig.12.These findings are summarized as the following observation.

    Table 2 Computation results of the DRO model.

    Fig. 9. Computation results of the DRO model. (a)α = 0.95; (b)α = 0.05.

    Fig. 10. Train schedule and in-vehicle passengers for the instance when α = 0.95, λ = 0.10, and Ψ = 0.02.

    Fig. 11. Convergence tendency of the best objective in the LS + MILP solver algorithm.

    Observation 2.(The necessity of considering λ depends on the risk attitude of decision-makers.) The numerical experiments in this section indicate the importance of the tradeoff coefficient λ in the proposed model, as it provides a support for generating a balanced solution that is customized to the actual needs of the decision-makers. The tradeoff coefficient λ provides decisionmakers with the opportunity to obtain a personalized and reliable plan.Decision-makers can choose the appropriate values of λ based on their subjective preferences and risk attitudes.

    5.3. Comparison between the SP model and the DRO model

    In this section,numerical experiments are conducted by solving the two-stage SP model (Eq. (27)). Since the probability distribution is known exactly in this model, the disturbance range of the uncertain probabilities—that is,Ψ—is equal to 0.Therefore,we display the computation results with different combinations of α and λ in Fig. 13. As shown, the optimal objective shows an increasing tendency with an increase in the risk-aversion parameter α if the tradeoff coefficient λ is fixed. In addition, when the level of risk aversion is fixed,the optimal objective follows the same trend with respect to the tradeoff coefficient λ.

    It should be noted that, since parameter Ψ is set as 0 in the SP model (Eq. (27)), the objective values should theoretically be less than those of the DRO model (Eq. (35)), with the same parameter settings. Below, we further analyze their differences between the SP model and the DRO model, by conducting a series of comparative experiments. For clarity, the following formula is used to quantify the differences(denoted as Price)in the optimal objective values between the SP model (termed SP*) and the DRO model(denoted as DRO*):

    Table 3 provides a detailed performance comparison for various combinations of α, λ, and Ψ.

    Fig. 12. Expected number of waiting passengers and CVaR with respect to λ. (a)α = 0.50,Ψ = 0.02; (b)α = 0.50,Ψ = 0.06.

    Fig. 13. Computation results of the SP model.

    From the computation results in Table 3,it can be observed that the DRO model always yields a slightly larger objective value than the SP model,since the values of Price are all positive.This demonstrates that, in an urban rail transit system, the more precise the information of probability distribution we know (i.e., the smaller Ψ is), the better the passenger flow control strategy and train schedule are. In this set of experiments, it can be seen that the maximum Price calculated is not greater than 0.051640%, which corresponds to a very small extra cost in the process of handling distribution uncertainty. Hence, in terms of solution stability, our proposed DRO model is superior to the SP model.

    Observation 3. (The superiority of the DRO model in comparison with the SP model.) To summarize, we conclude that the DRO model (Eq. (35)) utilizes limited distribution information to hedge against the uncertainty more robustly than the SP model(Eq. (27)), while requiring a relatively small Price.

    5.4. Comparison between the classic RO model and the DRO model

    In this section, we compare the performance of the classic RO model in Eq.(S2)and the DRO model(Eq.(35)).Since the risk measure, CVaR, is not involved in the classic robust method, we set λ = 0 in the following experiments. For concise descriptions, the classic RO model (Eq. (S2)) is referred to as the RO model. Fig. 14 provides a comparison of the results when adopting different Ψ and α in the models. As shown in Fig. 14, regardless of what Ψ is, the approximate optimal objective of the DRO model is smaller than that of the RO model. The reason for the poor optimization result of the RO model might be that the RO model is designed to optimize the worst-case scenario rather than the expected value, which leads to over-conservative solutions. It should be noted that Ψ represents the disturbance range of the uncertain probability.We thus conclude that the uncertainty can be handled well by our proposed DRO model,regardless of whether it is strong or weak.In addition,we set different values of α to carry out experiments, and the results show that the objective values remain the same when parameter α varies. This finding is consistent with the following analyses:Since α is set as zero in these experiments,the coefficient λ/(1-α)in the objective function is 0 for any value of α, so the corresponding part vanishes.

    Observation 4. (The superiority of the DRO model in comparison with the RO model.) From the computational results,we can derive the following observations: ①Our proposed DRO model outperforms the RO model in terms of solution qualityand stability. ②When a decision-maker is extremely risk averse,he/she can employ the classic RO method to optimize the worstcase scenario, so as to obtain a conservative solution; otherwise,we recommend the DRO approach with an appropriate combination of parameters, which can be selected flexibly according to the decision-makers’ preference and actual requirements.

    Table 3 The obtained Price in the numerical experiments.

    Fig. 14. Comparison of results from the RO model and the DRO model.

    5.5. Sensitivity analyses with respect to ζ1 and ζ2

    In this section, we perform sensitivity analyses with respect to the variation of two parameters, ζ1and ζ2. These two parameters are used in the objective function of the DRO model (Eq. (35)) to balance the operating time(denoted as Obj.1)related to the interests of the operators and the target relating to passengers(denoted as Obj.2).To be specific,if the value of ζ2/ζ1goes to ∞,then Obj.2 takes precedence over Obj. 1, which implies that more weight is given to passenger interests; otherwise, the proposed model mainly focuses on the optimization of the train schedule for operators.In particular,we aim to discuss how different values of ζ2/ζ1,leading to different train schedules,can affect the dynamic passenger control process.Table 4 reports the results obtained by varying ζ1,ζ2,and ζ2/ζ1.We list the objective values in the fifth column,and list the values of Obj.1 and Obj.2 in the sixth and seventh columns,respectively.As shown,when the value of ζ2/ζ1is between 0.01 and∞,the value of Obj.2 remains constant until ζ2/ζ1drops to 0.00001.The experiment with ζ1= 1 and ζ2= 0 does not consider the interests of the passengers, but only optimizes the train schedules. We note that the value of Obj. 2 increases drastically in this case. This computation result indicates the importance of collaboratively optimizing the robust passenger flow control strategies and the train schedule in order to ease congestion. In addition, in order to analyze the additional costs caused by the ambiguous distribution, we present the Price of each instance in Fig. 15. Clearly,although the value of Price varies considerably with respect to different ζ1and ζ2,it is still acceptable,since it does not exceed 4.00%at most. Moreover, for instances ζ1= 1 and ζ2= 0, the Price turns out to be 0, since the objective function does not take Obj. 2 into account.

    Fig. 15. The prices of distributional robustness under different values of ζ2/ζ1.

    Observation 5.(Sensitivity to parameters ζ1and ζ2.) Considering both Obj. 1 and Obj. 2, especially regarding the results of the last experiment, it can be seen that only optimizing the train schedule is insufficient to achieve the system optimum on the considered oversaturated line.To improve the service quality as much as possible, the decision-makers of the urban rail transit system must pay more attention to the collaborative optimization of robust passenger flow control and train scheduling approaches,and moreover, they must carefully determine the value of ζ2/ζ1based on the practical operation needs.

    6. Conclusions

    To reduce over-saturation of the urban rail transit line, this paper investigated the integrated problem of collaboratively optimizing passenger flow control strategies and train schedules, for a situation in which the passenger demand is uncertain and is characterized by a series of stochastic scenarios with partially known probability distribution. A multi-scenario two-stage DRO model was formulated to generate near-optimal strategies for the passenger flow control and train scheduling under distribution ambiguity. The mean-CVaR criterion was adopted in the objective function, and the goal was to minimize the sum of the operating time, the expected number of waiting passengers, and the CVaR of congestion, the coefficients of which could be adjusted flexibly according to the practical situation in order to achieve the most suitable balance. Since the proposed model is usually computationally intractable, a discrepancy-based ambiguity set was constructed by using partial information on the uncertain probability distribution of the stochastic scenarios. Computationally tractable forms of the proposed model were also deduced. Toobtain high-quality solutions for large-scale problems within an acceptable computing time, we developed a hybrid algorithm by combining an LS with a suitable MILP solver. To verify the performance of the proposed approaches, a series of numerical experiments were conducted using actual operation data from a Beijing rail transit line. The computation results show that our proposed DRO model is superior to the SP model in terms of solution stability. In addition, our model is superior in terms of solution quality and solution stability, compared with the classic RO model. These findings validate the proposed DRO method and the LS+MILP solver algorithm.

    Table 4 Sensitivity analyses of ζ1 and ζ2 with α = 0.50, λ = 0.10, and Ψ = 0.02.

    Future research can be dedicated to the following aspects:①In the process of characterizing the uncertain probability distribution,a discrepancy-based ambiguity set is used in this study. Thus, a promising extension would be to construct multi-type ambiguity sets for capturing the probability distribution, such as a Gaussian ambiguity set or a polyhedral ambiguity set. ②The study of integrated and interconnected multimodal transport systems remains a hot topic for the future development of urban mobile ecosystems.Hence,the collaborative optimization of the passenger flow control in a multimodal transportation system may be an interesting research direction for comprehensively easing congestion. ③In addition, it should be noted that this paper only considers the operation environment of an urban rail transit line. An interesting direction for future research would be to investigate efficient algorithms for more complex lines or networks in Beijing and in other representative cities.

    Acknowledgments

    This work was supported the National Natural Science Foundation of China (71621001, 71825004, and 72001019), the Fundamental Research Funds for Central Universities (2020JBM031 and 2021YJS203),and the Research Foundation of State Key Laboratory of Rail Traffic Control and Safety (RCS2020ZT001).

    Compliance with ethics guidelines

    Yahan Lu, Lixing Yang, Kai Yang, Ziyou Gao, Housheng Zhou,Fanting Meng, and Jianguo Qi declare that they have no conflict of interest or financial conflicts to disclose.

    Appendix A. Supplementary data

    Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2021.09.016.

    亚洲av成人一区二区三| 精品国内亚洲2022精品成人| 国产成人精品久久二区二区免费| 欧美性长视频在线观看| 黄色 视频免费看| 亚洲真实伦在线观看| 亚洲精品久久成人aⅴ小说| 婷婷精品国产亚洲av在线| 国产一区在线观看成人免费| 天天躁狠狠躁夜夜躁狠狠躁| 国产av一区二区精品久久| 成人亚洲精品一区在线观看| 男人的好看免费观看在线视频 | 国产精品一区二区三区四区久久 | 可以在线观看的亚洲视频| 成熟少妇高潮喷水视频| netflix在线观看网站| 777久久人妻少妇嫩草av网站| 伊人久久大香线蕉亚洲五| 一进一出好大好爽视频| 人成视频在线观看免费观看| 国产成人一区二区三区免费视频网站| 亚洲欧洲精品一区二区精品久久久| 国产熟女xx| 在线观看www视频免费| 手机成人av网站| 天天躁夜夜躁狠狠躁躁| 免费高清视频大片| 男人的好看免费观看在线视频 | 亚洲欧美精品综合久久99| 脱女人内裤的视频| 亚洲国产高清在线一区二区三 | 久久久久久国产a免费观看| 国产高清有码在线观看视频 | 亚洲,欧美精品.| 国产精品自产拍在线观看55亚洲| 51午夜福利影视在线观看| 男人的好看免费观看在线视频 | 法律面前人人平等表现在哪些方面| 成人精品一区二区免费| 91九色精品人成在线观看| 一区二区三区激情视频| 国产在线观看jvid| 免费看美女性在线毛片视频| 久久久精品国产亚洲av高清涩受| 黄片播放在线免费| 中文字幕久久专区| 久久精品91蜜桃| 91字幕亚洲| 色综合亚洲欧美另类图片| 国产一区二区三区在线臀色熟女| 午夜老司机福利片| 亚洲欧美日韩高清在线视频| 亚洲色图av天堂| 国产精品影院久久| 国产成年人精品一区二区| 久久久国产精品麻豆| 18禁美女被吸乳视频| 国产成人一区二区三区免费视频网站| 99精品在免费线老司机午夜| 日本熟妇午夜| avwww免费| 满18在线观看网站| 精品久久久久久久久久久久久 | www.精华液| 白带黄色成豆腐渣| 亚洲,欧美精品.| 嫁个100分男人电影在线观看| 国产久久久一区二区三区| 日本熟妇午夜| 男女下面进入的视频免费午夜 | 欧美最黄视频在线播放免费| 2021天堂中文幕一二区在线观 | 侵犯人妻中文字幕一二三四区| av欧美777| 一个人观看的视频www高清免费观看 | 免费无遮挡裸体视频| 男女下面进入的视频免费午夜 | 一本大道久久a久久精品| 村上凉子中文字幕在线| 叶爱在线成人免费视频播放| 午夜福利成人在线免费观看| 黄网站色视频无遮挡免费观看| 亚洲第一青青草原| 18禁黄网站禁片午夜丰满| 人人妻人人看人人澡| 国产精品久久久久久亚洲av鲁大| 成人午夜高清在线视频 | 亚洲国产欧洲综合997久久, | 97人妻精品一区二区三区麻豆 | 色播在线永久视频| 亚洲真实伦在线观看| 少妇裸体淫交视频免费看高清 | 久久久久久久久免费视频了| 日韩免费av在线播放| 亚洲中文av在线| 女警被强在线播放| 亚洲最大成人中文| 亚洲成人国产一区在线观看| 久99久视频精品免费| 久久国产亚洲av麻豆专区| 伊人久久大香线蕉亚洲五| 国产一区二区三区在线臀色熟女| √禁漫天堂资源中文www| 久热爱精品视频在线9| 宅男免费午夜| 99久久国产精品久久久| 两性午夜刺激爽爽歪歪视频在线观看 | 国产主播在线观看一区二区| 最近在线观看免费完整版| 成人欧美大片| 亚洲第一av免费看| 18美女黄网站色大片免费观看| 日韩欧美国产在线观看| www.自偷自拍.com| 亚洲一码二码三码区别大吗| 亚洲在线自拍视频| 无遮挡黄片免费观看| 国产高清有码在线观看视频 | 看片在线看免费视频| 嫩草影视91久久| 国内毛片毛片毛片毛片毛片| 亚洲第一电影网av| 精品国产超薄肉色丝袜足j| 老汉色av国产亚洲站长工具| 夜夜夜夜夜久久久久| 亚洲最大成人中文| 免费高清视频大片| 91麻豆av在线| 亚洲精品国产区一区二| 欧美激情高清一区二区三区| 伦理电影免费视频| 999精品在线视频| 日韩三级视频一区二区三区| 少妇的丰满在线观看| 婷婷亚洲欧美| 波多野结衣高清作品| 日韩高清综合在线| 黑人欧美特级aaaaaa片| 亚洲精品一区av在线观看| 国产黄色小视频在线观看| 美女扒开内裤让男人捅视频| 999久久久国产精品视频| 听说在线观看完整版免费高清| 精品日产1卡2卡| 亚洲欧美精品综合一区二区三区| 视频区欧美日本亚洲| 国产精品亚洲av一区麻豆| 激情在线观看视频在线高清| 哪里可以看免费的av片| 亚洲精品国产一区二区精华液| 黑人巨大精品欧美一区二区mp4| 亚洲久久久国产精品| 久久这里只有精品19| 熟妇人妻久久中文字幕3abv| 国产精品二区激情视频| 老司机午夜福利在线观看视频| 亚洲av五月六月丁香网| av有码第一页| 日本三级黄在线观看| 亚洲 欧美 日韩 在线 免费| 免费看美女性在线毛片视频| 欧美一级a爱片免费观看看 | 亚洲国产精品成人综合色| 91老司机精品| 久久久国产欧美日韩av| 日本 欧美在线| 国产激情欧美一区二区| 亚洲片人在线观看| 真人一进一出gif抽搐免费| 亚洲av成人一区二区三| 亚洲精品色激情综合| 美女 人体艺术 gogo| 欧美精品啪啪一区二区三区| 天堂√8在线中文| 成年免费大片在线观看| 久久久久久免费高清国产稀缺| 90打野战视频偷拍视频| 日韩欧美在线二视频| 久久久久久人人人人人| 久久久久久久精品吃奶| 中文字幕高清在线视频| 欧美午夜高清在线| 在线av久久热| 久久婷婷成人综合色麻豆| 叶爱在线成人免费视频播放| 50天的宝宝边吃奶边哭怎么回事| 国产一区在线观看成人免费| 国产国语露脸激情在线看| 中文字幕精品免费在线观看视频| 国产av不卡久久| 日本撒尿小便嘘嘘汇集6| 白带黄色成豆腐渣| 欧美日韩中文字幕国产精品一区二区三区| 欧洲精品卡2卡3卡4卡5卡区| 国产av在哪里看| 久久精品91无色码中文字幕| 亚洲专区字幕在线| 国产av不卡久久| 亚洲第一欧美日韩一区二区三区| 老司机靠b影院| 大型av网站在线播放| 国产又黄又爽又无遮挡在线| 亚洲国产精品999在线| 99riav亚洲国产免费| 一区二区三区国产精品乱码| 国产精品,欧美在线| 国产极品粉嫩免费观看在线| 精品不卡国产一区二区三区| 精品免费久久久久久久清纯| 国产99久久九九免费精品| 亚洲欧美日韩无卡精品| av福利片在线| 国产亚洲精品久久久久久毛片| 久久久国产精品麻豆| 欧洲精品卡2卡3卡4卡5卡区| 天天躁夜夜躁狠狠躁躁| 高清毛片免费观看视频网站| 欧美性长视频在线观看| 亚洲狠狠婷婷综合久久图片| 精品国产乱子伦一区二区三区| 丁香欧美五月| 午夜精品在线福利| 国产精品98久久久久久宅男小说| 国产亚洲精品久久久久5区| 久久久久久久久免费视频了| 草草在线视频免费看| 两性夫妻黄色片| 欧美国产精品va在线观看不卡| 国产精品二区激情视频| 91麻豆av在线| 国产三级在线视频| 久99久视频精品免费| 欧美精品亚洲一区二区| 久久久久久国产a免费观看| 午夜成年电影在线免费观看| 亚洲国产精品sss在线观看| 黄片大片在线免费观看| 老司机在亚洲福利影院| 757午夜福利合集在线观看| 亚洲熟妇中文字幕五十中出| 日本三级黄在线观看| 午夜福利在线在线| 在线观看日韩欧美| 黄色成人免费大全| 日本a在线网址| 日韩大尺度精品在线看网址| 1024香蕉在线观看| 最新在线观看一区二区三区| 一进一出抽搐gif免费好疼| 69av精品久久久久久| 真人做人爱边吃奶动态| 岛国视频午夜一区免费看| 国产v大片淫在线免费观看| 老司机靠b影院| 亚洲精华国产精华精| 亚洲av成人一区二区三| 美女高潮到喷水免费观看| 精品久久久久久久久久久久久 | 在线观看日韩欧美| 久久久久久免费高清国产稀缺| 欧美zozozo另类| 波多野结衣巨乳人妻| 宅男免费午夜| 欧美中文综合在线视频| 国产在线观看jvid| 亚洲一区高清亚洲精品| 国产精品久久久久久亚洲av鲁大| 成人手机av| 欧美乱码精品一区二区三区| АⅤ资源中文在线天堂| 色尼玛亚洲综合影院| 自线自在国产av| 国产精品久久久久久人妻精品电影| 精品国产乱子伦一区二区三区| 午夜福利欧美成人| 午夜精品久久久久久毛片777| 成人国产综合亚洲| 听说在线观看完整版免费高清| 久久热在线av| 一a级毛片在线观看| 国产精华一区二区三区| www日本在线高清视频| 亚洲,欧美精品.| 亚洲成a人片在线一区二区| 女同久久另类99精品国产91| 99国产精品一区二区蜜桃av| 淫妇啪啪啪对白视频| 天天躁狠狠躁夜夜躁狠狠躁| 后天国语完整版免费观看| 搡老熟女国产l中国老女人| 一进一出好大好爽视频| 欧美成人免费av一区二区三区| 亚洲精品国产区一区二| 在线观看午夜福利视频| 日本熟妇午夜| 天天添夜夜摸| 精品国产乱子伦一区二区三区| 国产一区在线观看成人免费| 97超级碰碰碰精品色视频在线观看| 亚洲专区国产一区二区| 两个人视频免费观看高清| 亚洲欧美日韩高清在线视频| 他把我摸到了高潮在线观看| 成熟少妇高潮喷水视频| 亚洲国产精品sss在线观看| 热99re8久久精品国产| 国产一级毛片七仙女欲春2 | 日本五十路高清| 后天国语完整版免费观看| 少妇 在线观看| а√天堂www在线а√下载| 精品国产国语对白av| 亚洲av成人av| 精品国内亚洲2022精品成人| АⅤ资源中文在线天堂| 精品第一国产精品| 国产91精品成人一区二区三区| 99热这里只有精品一区 | 女性被躁到高潮视频| 99国产综合亚洲精品| 丁香六月欧美| 色综合站精品国产| 人人妻人人看人人澡| 欧美激情极品国产一区二区三区| 国产99久久九九免费精品| 精品乱码久久久久久99久播| www.熟女人妻精品国产| 日韩 欧美 亚洲 中文字幕| 首页视频小说图片口味搜索| 色在线成人网| 可以免费在线观看a视频的电影网站| 久久久久亚洲av毛片大全| 真人一进一出gif抽搐免费| 久久香蕉精品热| а√天堂www在线а√下载| 国产伦一二天堂av在线观看| 欧美激情 高清一区二区三区| 国产精品永久免费网站| 精品国产美女av久久久久小说| 99久久综合精品五月天人人| 国产又爽黄色视频| 99久久综合精品五月天人人| 亚洲国产看品久久| 国产亚洲欧美精品永久| 亚洲国产看品久久| 国产久久久一区二区三区| 国产又色又爽无遮挡免费看| 午夜久久久在线观看| 啦啦啦免费观看视频1| 狂野欧美激情性xxxx| 午夜福利18| 国产成人一区二区三区免费视频网站| 欧美又色又爽又黄视频| 成年免费大片在线观看| 丁香欧美五月| 国产激情久久老熟女| 搡老岳熟女国产| 在线av久久热| 精品国产一区二区三区四区第35| 亚洲国产中文字幕在线视频| 久热这里只有精品99| 少妇的丰满在线观看| 无遮挡黄片免费观看| 欧美黄色淫秽网站| 亚洲第一av免费看| 法律面前人人平等表现在哪些方面| av在线播放免费不卡| 午夜福利成人在线免费观看| 国产av一区二区精品久久| 国产久久久一区二区三区| 国产高清激情床上av| 欧洲精品卡2卡3卡4卡5卡区| 亚洲激情在线av| 国产精品永久免费网站| 一进一出好大好爽视频| 亚洲在线自拍视频| 香蕉丝袜av| 精品国产国语对白av| 听说在线观看完整版免费高清| 91在线观看av| 日本a在线网址| 又黄又爽又免费观看的视频| 久久天躁狠狠躁夜夜2o2o| 一二三四在线观看免费中文在| 日本 欧美在线| 一级黄色大片毛片| 欧美乱妇无乱码| 国产一区二区三区视频了| 欧美色视频一区免费| 亚洲成人免费电影在线观看| 欧美人与性动交α欧美精品济南到| 老汉色∧v一级毛片| 午夜精品久久久久久毛片777| 国产爱豆传媒在线观看 | 色综合婷婷激情| 国产精品自产拍在线观看55亚洲| 99国产精品一区二区三区| 日日干狠狠操夜夜爽| 欧美精品啪啪一区二区三区| 久久天躁狠狠躁夜夜2o2o| 国产亚洲av高清不卡| 91国产中文字幕| 免费在线观看影片大全网站| 免费在线观看亚洲国产| 在线视频色国产色| 欧美黄色淫秽网站| 日本 av在线| 久久久久精品国产欧美久久久| 国产精品,欧美在线| 欧美一级a爱片免费观看看 | e午夜精品久久久久久久| 亚洲欧美日韩无卡精品| 日本三级黄在线观看| 亚洲激情在线av| 黄色a级毛片大全视频| 国产一区二区在线av高清观看| 亚洲中文日韩欧美视频| 免费看美女性在线毛片视频| 亚洲色图av天堂| 日韩高清综合在线| www日本黄色视频网| 亚洲第一欧美日韩一区二区三区| 黄色女人牲交| 亚洲在线自拍视频| 成人亚洲精品一区在线观看| 两性午夜刺激爽爽歪歪视频在线观看 | 夜夜爽天天搞| 欧美三级亚洲精品| 亚洲人成电影免费在线| 村上凉子中文字幕在线| 日本三级黄在线观看| 亚洲激情在线av| 精品一区二区三区av网在线观看| 淫秽高清视频在线观看| 欧美精品啪啪一区二区三区| 18禁裸乳无遮挡免费网站照片 | 免费搜索国产男女视频| 亚洲成人免费电影在线观看| 久9热在线精品视频| 99久久久亚洲精品蜜臀av| 午夜精品久久久久久毛片777| 老司机靠b影院| 很黄的视频免费| 日韩成人在线观看一区二区三区| 丰满人妻熟妇乱又伦精品不卡| 国产91精品成人一区二区三区| 欧美久久黑人一区二区| 成人亚洲精品一区在线观看| 色婷婷久久久亚洲欧美| 丝袜人妻中文字幕| 不卡一级毛片| tocl精华| av天堂在线播放| 国产av一区在线观看免费| 午夜免费观看网址| 99国产精品99久久久久| 亚洲av第一区精品v没综合| 欧美精品亚洲一区二区| 国产真实乱freesex| 男女午夜视频在线观看| 欧美国产日韩亚洲一区| 亚洲熟妇熟女久久| 香蕉久久夜色| 欧美黑人巨大hd| 欧美日韩黄片免| 无人区码免费观看不卡| av在线播放免费不卡| 久久国产精品人妻蜜桃| or卡值多少钱| 757午夜福利合集在线观看| 久久久久九九精品影院| 欧美又色又爽又黄视频| www.精华液| 十八禁网站免费在线| 中文资源天堂在线| 老鸭窝网址在线观看| 国产成人欧美| 欧美黄色淫秽网站| 91麻豆精品激情在线观看国产| 黄色毛片三级朝国网站| 亚洲七黄色美女视频| 欧美不卡视频在线免费观看 | 精品国产乱子伦一区二区三区| 亚洲狠狠婷婷综合久久图片| 国产精品乱码一区二三区的特点| 久久精品国产亚洲av高清一级| 又黄又粗又硬又大视频| 日本免费a在线| 精品国产国语对白av| 美女高潮喷水抽搐中文字幕| 亚洲aⅴ乱码一区二区在线播放 | 久久久久国内视频| av免费在线观看网站| 亚洲七黄色美女视频| 亚洲午夜精品一区,二区,三区| 满18在线观看网站| 国产乱人伦免费视频| 国产精品 国内视频| 国产精品1区2区在线观看.| 色播亚洲综合网| 色婷婷久久久亚洲欧美| 人成视频在线观看免费观看| 手机成人av网站| 母亲3免费完整高清在线观看| 色老头精品视频在线观看| 亚洲黑人精品在线| 欧美不卡视频在线免费观看 | 桃色一区二区三区在线观看| 久久人人精品亚洲av| 午夜激情福利司机影院| 婷婷精品国产亚洲av在线| 亚洲欧美激情综合另类| 此物有八面人人有两片| 欧美色视频一区免费| 国产精品综合久久久久久久免费| 中文字幕人妻熟女乱码| 久久久久久九九精品二区国产 | 人妻久久中文字幕网| 精华霜和精华液先用哪个| 午夜福利成人在线免费观看| 日韩 欧美 亚洲 中文字幕| 欧美另类亚洲清纯唯美| 午夜a级毛片| 黄色丝袜av网址大全| 久久天堂一区二区三区四区| 在线播放国产精品三级| 一进一出好大好爽视频| 人人澡人人妻人| 十八禁网站免费在线| 久久 成人 亚洲| 久久久精品欧美日韩精品| 欧美乱码精品一区二区三区| 国产伦人伦偷精品视频| 日韩精品中文字幕看吧| 久久精品91无色码中文字幕| 黄网站色视频无遮挡免费观看| 桃红色精品国产亚洲av| 男人舔奶头视频| 黄色丝袜av网址大全| 日韩视频一区二区在线观看| 侵犯人妻中文字幕一二三四区| 亚洲人成伊人成综合网2020| 午夜免费观看网址| 国产成+人综合+亚洲专区| 18禁裸乳无遮挡免费网站照片 | 免费在线观看黄色视频的| av欧美777| 精品高清国产在线一区| av超薄肉色丝袜交足视频| 日本免费一区二区三区高清不卡| 亚洲国产精品久久男人天堂| 国产一区二区三区在线臀色熟女| 午夜两性在线视频| 久久精品国产99精品国产亚洲性色| 在线免费观看的www视频| 中文资源天堂在线| 禁无遮挡网站| 国产一区在线观看成人免费| 黄色毛片三级朝国网站| 18美女黄网站色大片免费观看| av天堂在线播放| 两个人看的免费小视频| 少妇裸体淫交视频免费看高清 | 免费在线观看完整版高清| 中文资源天堂在线| 欧美乱妇无乱码| 色在线成人网| 精品久久久久久久久久免费视频| 亚洲七黄色美女视频| ponron亚洲| 久久草成人影院| 亚洲精品美女久久av网站| 久久久久国产精品人妻aⅴ院| 麻豆成人午夜福利视频| 婷婷六月久久综合丁香| 人人妻人人看人人澡| 一区二区三区精品91| 亚洲国产毛片av蜜桃av| 黄片播放在线免费| 亚洲aⅴ乱码一区二区在线播放 | videosex国产| 国产视频一区二区在线看| 久久精品国产99精品国产亚洲性色| 久久久久久久久久黄片| 国产亚洲精品av在线| 高清毛片免费观看视频网站| 精华霜和精华液先用哪个| 麻豆一二三区av精品| 巨乳人妻的诱惑在线观看| 脱女人内裤的视频| 中亚洲国语对白在线视频| 久久欧美精品欧美久久欧美| 最近最新中文字幕大全电影3 | 国产精品野战在线观看| 亚洲精品久久国产高清桃花| 又紧又爽又黄一区二区| 欧美中文日本在线观看视频| 久久精品影院6| 国产精品免费视频内射| 性欧美人与动物交配| 天堂动漫精品| av免费在线观看网站| 十分钟在线观看高清视频www| 天堂√8在线中文| 国产精品一区二区精品视频观看| 国产97色在线日韩免费| 黄色视频不卡| 国产亚洲精品久久久久久毛片| 欧美成人免费av一区二区三区| 露出奶头的视频| 999精品在线视频| 国产av一区在线观看免费| 国产亚洲精品久久久久5区| 999久久久国产精品视频|