Modeling and Identification of Nonlinear Dynamics for Freeway Traffic ...
952
IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 17, NO. 4, JULY 2009
Modeling and Identification of Nonlinear Dynamics for Freeway Traffic by Using Information From a Mobile Cellular Network Angelo Alessandri, Senior Member, IEEE, Raffaele Bolla, Member, IEEE, Mauro Gaggero, and Matteo Repetto
Abstract—The high coverage of the territory by cellular networks and the widespread diffusion of mobile terminals aboard vehicles allow one to collect information on the traffic behavior. The problem of selecting a dynamic model to describe the freeway traffic by using the information available from a wireless cellular network is addressed by assuming the distribution of mobile terminals aboard vehicles to be uniform along the carriageway. Two different nonlinear parametrized models of freeway traffic are investigated: the first is an extension to a well-established macroscopic model, while the second is based on a black-box approach and consists in using a neural network to approximate the traffic dynamics. The parameters of such models are identified off line by a least-squares technique. Traffic measurements obtained from a cellular network are employed to identify and validate the proposed models, as shown by means of simulations. Index Terms—Freeway traffic model, identification, least squares, macroscopic model, mobile cellular network, neural networks.
I. INTRODUCTION N THIS brief, Dynamic models to describe the traffic on a freeway stretch are proposed that are based on the information about density, mean speed, and flow of vehicles that can be provided by a wireless cellular network, with a limited use of traditional measurement devices such as loop detectors, radars, and television cameras (see, e.g., [1]). Since such sensors are very expensive, a reduced number of them is placed along the roads, usually at long distances, thus providing only local information, possibly affected by inhomogeneities [2]. Here we shall exploit the traffic information that can be obtained from mobile terminals aboard vehicles, i.e., phones, handhelds, and personal computers in ongoing communications with a base station for a voice call or a data transfer (we shall denote such terminals as “active”). Starting from the ideas reported in [3], a section between two standard sensors is replaced with a cell of a wireless network. The cellular network technologies allow one to have a large number of small cells, 2–4 km in diameter, as compared with much larger distances between successive sensors. Moreover, the diffusion of wireless communications makes it very likely for vehicles to take aboard mobile terminals. The main
I
Manuscript received July 11, 2008; revised November 13, 2008. Manuscript received in final form January 22, 2009. First published May 19, 2009; current version published June 24, 2009. Recommended by Associate Editor C. C. Lee. A. Alessandri and M. Gaggero are with the Department of Production Engineering, Thermoenergetics, and Mathematical Models (DIPTEM), University of Genoa, I-16129 Genoa, Italy (e-mail: [email protected]; [email protected]). R. Bolla and M. Repetto are with the Department of Communications, Computer and System Sciences (DIST), University of Genoa, I-16145 Genoa, Italy (e-mail: [email protected]; [email protected]). Digital Object Identifier 10.1109/TCST.2009.2014242
drawbacks of standard acquisition devices (i.e., locality of traffic measurements and high costs) are overcome by acquiring information from mobile devices. Using the cellular network, one can directly collect the measurements of density of vehicles with active terminals aboard and flow of such vehicles moving from a cell to the next one. Moreover, recently developed techniques (see, e.g., [4]) allow one to determine the speed of active terminals, and hence of vehicles having them aboard. It is worth noting that, if required, one can integrate the information obtained from a cellular network with those available from traditional traffic sensors. The idea of collecting traffic measurements from mobile terminals aboard vehicles is reported also in [5], where the problem of estimating density and mean speed of vehicles is addressed by using a dynamic description of the traffic based on a partial differential equation that accounts for the conservation of the density of vehicles. In this respect, our contribution is that of developing two parametrized models that include also the speed dynamics and focusing on the identification of their parameters. Under the assumption on the uniform distribution of mobile terminals inside each cell of a wireless network, we propose two different approaches to model the freeway traffic. The first is an extension of a well-known macroscopic traffic model (see, e.g., [2], [6]–[8]). Such extension consists in adding a state variable that represents the percentage of vehicles with on-board active terminals with respect to the total number of vehicles in each cell. The second relies on the idea of using a neural network to approximate the traffic dynamics according to the black-box paradigm, and hence does not require all the fluid hypotheses used to construct the first. Such an approach is motivated by the excellent approximating properties of neural networks and the availability of efficient methods to tune the neural weights (see, e.g., [9]–[13]). Both the first and the second model are nonlinear in the state variables and depend nonlinearly on the corresponding parameters, which can be identified using least squares. In order to validate the models that result from the identification of the corresponding parameters, a simulation setup has been devised that is composed of: 1) a freeway traffic microscopic simulator that traces the behaviors of all the vehicles moving on the freeway (see [14], [15]) and 2) a cellular network simulator that takes into account voice calls or data transfers generated by on-board terminals. The rest of the brief is organized as follows. Section II describes the proposed extended macroscopic traffic model and the least-squares approach used for the identification of its parameters. In Section III, the same problems of modeling and identification are addressed via neural networks. Section IV reports the results obtained in a simulation case study. The conclusions are drawn in Section V.
1063-6536/$25.00 © 2009 IEEE
ALESSANDRI et al.: MODELING AND IDENTIFICATION OF NONLINEAR DYNAMICS FOR FREEWAY TRAFFIC
953
Fig. 1. Macroscopic model of traffic in a freeway trunk.
II. EXTENDED MACROSCOPIC TRAFFIC MODEL In the literature, there exist two main families of models to represent the traffic behavior along a road: the macroscopic ones, which deal only with aggregate quantities (e.g., density and mean speed of vehicles), and the microscopic ones, which describe the behavior of all the single vehicles. Using a microscopic model, one can easily simulate particular events such as slowing-downs, accidents, and carriageway restrictions. A microscopic simulator may be very detailed and precise, but requires a great computational effort. By contrast, in a macroscopic framework, the traffic is modeled like a fluid by using continuous variables, thus turning in less flexibility for simulating specific situations like the aforementioned ones; anyway, as compared with the microscopic approach, the computational burden is reduced. A macroscopic model of the traffic on a freeway trunk can be obtained by dividing the stretch into various sections, each described by state variables such as density and mean speed of vehicles (see Fig. 1). Moreover, one has to take into account the traffic flows leaving a section and entering the following one, as well as the input and output flows at the on- and offramps, respectively. We assume to perform a discretization at time instants , where and is the sampling period. At time , the following variables are defined in a freeway stretch made up of segments: 1) is the density of vehicles in section (in veh/km), ; 2) is the mean speed of vehicles in section (in km/h), ; 3) is the flow of vehicles leaving section and entering section (in veh/h), ; 4) and are the flows of vehicles entering and coming out of section , by means of on- and off-ramps (in veh/h), respectively; and are measured via standard sensors located at the ramps and are identically null if there exist no ramps in section . Second-order models are usually considered to describe the dynamics of freeway traffic by using such aggregate variables. Here we shall adopt the Payne–Whitham model [16], [17], by referring to the well-established approach developed by Papageorgiou and coworkers (see [2], [6]–[8], and the references therein). The density of vehicles in section varies from a time interval to the next, depending on the flows of vehicles that entered and left the section, according to the following balance equation:
(1)
, and is the length of section (in km), where which is assumed to be known. For the first and last sections, state-augmentation equations are considered, i.e., and , where the scalar noises and model the uncertainties in the densities at the boundaries. Concerning the speed dynamics, the following nonlinear equation is adopted:
(2) where and is a noise. The speed dynamics in the boundary segments is described by and , where and are scalar noises. Moreover, , and are positive constants and is the function representing the static speed-density relationship, where is the free speed, and is the critical density. Finally, the flow of vehicles between section and section is given by (3) and is a constant parameter. where The model given by (1), (2), and (3) can be extended by using the traffic measurements from a cellular network that refer to a subset of vehicles (i.e., those with an active terminal aboard). Such information is basically the density of active terminals in each cell of the network and the handoff flows. Moreover, methodologies are now available to estimate the speed of a moving active terminal by using various techniques (see [4] and the references therein). Thus, we shall present a model that allows one to predict the behavior of the overall traffic in a freeway stretch by using the information obtained from a cellular network. Such information refers to a subset of vehicles (i.e., density, mean speed, and flow only of the vehicles with an active terminal aboard). If just one mobile terminal were on each vehicle and one would know the total number of such terminals in each cell of the wireless network, the density of vehicles on each cell could be computed exactly. Unfortunately, the cellular networks do not know the precise cell inside which each terminal resides. As a matter of fact, in order to limit the signalling between the
954
IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 17, NO. 4, JULY 2009
Fig. 2. Extended macroscopic model of traffic in a freeway trunk.
mobile terminals and the network, the positions of the terminals are updated only with respect to the location area, which is a set of quite a large number of cells that represents the most precise information about the positions of the terminals. However, when a mobile terminal is active (i.e., when it is communicating with a base station), the cell on which it resides is known. In other words, a cellular network can supply information about the number of active terminals on each cell and the flows between adjacent cells due to handoffs that occur because of the movement of active terminals from one cell to the next one along the freeway. In the following, we denote the vehicles with an active terminal aboard as the “observed” ones. The freeway sections are replaced with network cells. The ratio between the number of vehicles with active terminals aboard and the total number of vehicles in cell at time is denoted by and regarded as a state variable for each cell together with density and mean speed. The resulting model is depicted in Fig. 2. To sum up, the description of the proposed model requires the definition of the following additional variables at time : 1) is the percentage of observed vehicles with respect to the total number of vehicles in cell ; 2) is the density of observed vehicles in cell (in veh/ km), ; 3) is the density variation of observed vehicles in cell (in veh/km) due to the entering/leaving from the ramps or the start/end of voice calls or data transfers, ; 4) is the flow of observed vehicles leaving cell and entering cell (in veh/h), . The balance of observed vehicles entering and leaving cell yields the following relationship:
where
, and (6)
, and is given by (3). Thus, where can be computed by means of (4) and (6) as follows:
(7) where . For the boundary segments, the percentages of observed vehicles are assumed to satisfy the following relationships: and , where and are scalar noises. To derive a dynamic equation for that accounts for the measurements provided by the cellular network, we substitute (6) into (1) so as to obtain
(8) where . We recall that the measures of the variables , and are available from the cellular network, as well as and can be easily determined by standard sensors located at the on- and off-ramps, respectively. Moreover, we assume to measure the mean speed of vehicles by averaging the speeds of active terminals; we denote these measures by . Likewise in [5], we suppose that the mean speeds of the observed vehicles and of all the vehicles approximately agree under the uniformity assumption, i.e., (9)
where dynamics for
. Using this equation and (1), we derive the , i.e.,
(4) where . If we assume the distribution of vehicles with active terminals to be uniform inside each cell, the following relationships hold: (5)
(note that in the following we shall introduce where measurement noises to model uncertainties in the measures). Therefore, keeping in mind (5), (6), and (9), let us define the measurement vector as
where
and, from now on, the notation stands for a generic column vector . In addition, let us introduce the following definitions at time : 1) is the vector of the state variables;
ALESSANDRI et al.: MODELING AND IDENTIFICATION OF NONLINEAR DYNAMICS FOR FREEWAY TRAFFIC
2) is the vector that collects all the uncertainties affecting the dynamics; 3) is the vector of the density variations of observed vehicles due to their entering/leaving from the ramps or the start/end of voice calls or data transfers; 4) and are the vectors containing the flows of vehicles at the on- and off-ramps, respectively. For the sake of compactness, let us denote the difference equation and the output mapping that represent the dynamic equations [i.e., (8), (2), and (7)] and measurement equations [i.e., (5), (6), and (9)] by means of the functions and , respectively: (10a) (10b) where is the vector of the parameters to be determined together with , and is the vector of the measurement noises. The model (10) depends nonlinearly on the parameters and . According to a classical least-squares approach (see, e.g., [18]), the noise vectors and are supposed to be “small,” so one can try to estimate and by minimizing a regression cost function based on a batch of data given by the measures , and for . On the basis of the foregoing, after choosing , the predicted state is given by
955
rameters (i.e., neural weights) to ensure a fixed approximation accuracy, especially in high-dimensional settings. More specifically, one-hidden-layer sigmoidal neural networks and radialbasis-function networks with tunable external and internal parameters enjoy the property of uniform approximation and guarantee a precision that grows at most polynomially with the dimension of the input of the function to be approximated (see, e.g., [9]–[12] and the references therein). It is worth noting that such a property does not hold for linear approximators. Another crucial aspect that may suggest the use of neural approximators is the possibility of using efficient methodologies for optimizing the neural weights (e.g., the Levenberg–Marquardt or Newton algorithms [13]). For the purpose of comparison with the extended macroscopic model proposed in Section II, we consider a model of the system (10) without all the fluid assumptions introduced to derive the function . Thus, an alternative expression of is searched for by using a neural mapping and the available measurements given by . Clearly, this entails the necessity of relying on the assumptions used to construct the function . Toward this end, according to well-established approaches available in the literature (see, among others, [19]–[21]), let us consider the black-box parametrized model
where
and the function
is a neural mapping that is parametrized by the -dimensional neural weights vector with the compact set . As in Section II, let us introduce a least-squares cost function that depends on and , i.e.,
(11) where and is some prediction of the initial state. A least-squares quadratic cost function that measures the quality of the fitting is
(12) where is a symmetric, positive definite matrix. An estimate of and can be obtained by minimizing the cost (12) under the constraints regarding the dynamic equations and the ranges of values that can be taken on by the parameters and the state variables. Thus, the problem reduces to search for and that minimize (12) under the constraints (11), . III. BLACK-BOX MODELING USING NEURAL NETWORKS The use of neural networks for the purpose of modeling is motivated by the fact that continuous functions can be approximated arbitrarily well on a compact set by large classes of neural networks with one or more hidden layer(s). Such a capability is referred to as “universal approximation property.” Besides, one-hidden-layer neural networks exhibit another powerful feature that consists in requiring a small number of pa-
(13) where is a symmetric, positive definite matrix. The optimal weights and parameter can be estimated by minimizing the cost (13) under the same type of constraints considered in the minimization of the cost (12) for the identification of the extended macroscopic model. IV. SIMULATION RESULTS A. Description of the Simulation Setup In order to evaluate the effectiveness of the proposed models, a simulation setup has been devised that is made up of a microscopic freeway traffic simulator and a cellular network simulator. The microscopic traffic simulator is based on the car-following model proposed by Gipps and adapted to the freeway environment [14], [15]. The Gipps model allows one to describe the typical driver behavior according to the physical laws of motion by taking into account the traffic conditions around each vehicle. The main idea behind such a modeling approach is that each driver reacts to a stimulus coming from other vehicles on the road, especially the previous and the following ones. Different situations of possible interest can be easily simulated, like, for example, accidents and lane restrictions, which are
956
IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 17, NO. 4, JULY 2009
Fig. 3. Estimates of the parameters of the extended macroscopic model during the optimization process.
Fig. 4. Behaviors of the state variables of the extended macroscopic model compared with the corresponding ones given by the microscopic traffic and cellular network simulators.
particularly well-suited to testing our models. We used the microscopic traffic simulator to reproduce the real behavior of vehicles on a freeway. The data generated by such simulator together with the cellular network simulator are used as measurements to compute the state variables of the macroscopic and neural models. The values of the percentages of observed vehicles, density, and mean speed of all vehicles given by such simulators are used for comparison with the same variables obtained by the proposed models.
As to the cellular network simulator, we considered the typical operations of a generic wireless network at the connection level since no details concerning the protocol (e.g., ETACS, GSM, UMTS, WiFi, WiMax) and the switching technique (circuit or packet) are taken into account. Radio resources are limited for each cell of the network, and hence the number of concurrent connections is limited as well. In the following, we shall focus on a simulation case study consisting of a one-way freeway stretch of 24 km covered by
ALESSANDRI et al.: MODELING AND IDENTIFICATION OF NONLINEAR DYNAMICS FOR FREEWAY TRAFFIC
Fig. 5. Autocorrelation functions of the residuals of q
957
(t); i = 4; 5; 6, obtained with the extended macroscopic model.
Fig. 6. Behaviors of some state variables of the neural black-box model compared with the corresponding ones given by the microscopic traffic and cellular network simulators.
12 cells of a cellular network (numbered sequentially from 0 to 11). Each cell is 2 km long, and there exist on- and off-ramps in cells 2, 6, and 10. The reference speed of each vehicle, the rate of entrance from the beginning of the road, as well as the rates of entrance and exit from the ramps, were randomly chosen according to uniform distributions in the ranges [70, 100] km/h, [500, 1000] veh/h, [230, 270] veh/h, and [170, 270] veh/h, respectively. The connections to the network were generated according to a Poisson distribution with mean proportional to the vehicle density [the proportional coefficient was chosen to be equal to 150 Erlang km/(s veh)]. The sampling period was taken equal to 10 s and each simulation was 2 h long.
For a quantitative evaluation of the model fitting, we applied standard autocorrelation tests on the residuals between the measured variables and their predictions given by the outputs of a model (see [22] for an introduction). Toward this end, let us define the residuals of a scalar measurement sequence , as the difference between and its corresponding prediction . The autocorrelation function of the residuals , is defined as
958
IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 17, NO. 4, JULY 2009
Fig. 7. Autocorrelation functions of the residuals of q
where
and .
(t); i = 4; 5; 6, obtained with the neural black-box model.
is the sample mean, i.e.,
B. Identification and Validation of the Extended Macroscopic Model In order to identify the parameters of the extended macroscopic model, we used the routine fmincon of the Matlab Optimization Toolbox (see [23] for details) to minimize the cost (12) with equal to and a batch of data that takes into consideration various events occurring on a freeway stretch (e.g., slowing-downs and accidents). The estimates of the parameters at each iteration of the optimization are shown in Fig. 3. The parameters converged asymptotically to the following values: 0.0031 h, 40.56 km /h, 47.97 veh/km, 105.14 km/h, 70.48 veh/km, and . For the sake of brevity, we report only the simulation results with the behaviors of the state variables of the identified model in the central cells (i.e., cells 4–6), with a slowing-down occurring in cell 4 (see Fig. 4). The dashed lines represent the “real-world values” of the state variables given by the microscopic traffic and cellular network simulators, while the continuous lines represent the state variables of the extended macroscopic model with the identified parameters. In Fig. 5, the autocorrelation functions of the residuals of , are shown, where the horizontal lines bound the 95% confidence bands. The autocorrelation functions lie outside the 95% confidence bands for a small number of lags, thus showing a satisfying behavior of the identified model. C. Identification and Validation of the Neural Black-Box Model We used one-hidden-layer feedforward neural networks with different numbers of sigmoidal units in the hidden layer to identify only the dynamics of the central cells (i.e., cells 4–6). Therefore, the neural networks had only 21 inputs and 9 outputs, which considerably reduced the computational burden for the selection of the weights and parameter . The minimization of the cost (13) with equal to was performed using the routine fmincon of the Matlab Optimization Toolbox. Among the various attempts, Fig. 6 shows the best results obtained by a one-hidden-layer feedforward neural network with 10 sigmoidal units in the hidden layer. The value of after the identification turned out to be equal to 0.8. By comparing such results with those obtained by the extended macroscopic
model (see Fig. 4), one can see that the speed dynamics is identified very well, while the percentage of on-board active terminals suffer from the problem of a small bias. As to the density of vehicles, the performances of the two models are quite similar. However, the autocorrelation functions of the residuals of , given by the neural black-box model lie outside the 95% confidence bands for a larger number of lags with respect to those of the extended macroscopic model, as shown in Fig. 7. V. CONCLUSION The problem of modeling and identifying the dynamics of freeway traffic by using the information available from active mobile terminals aboard vehicles has been faced according to two different paradigms. The first results from an extension to a well-known macroscopic model of freeway traffic. The second is based on the idea of using a black-box neural model, thus avoiding a lot of fluid assumptions that are required by the first approach. The analysis of the results obtained by the identification shows that the former approach has turned out to perform better than the latter on the overall. Future research will concern the use of the proposed models for the purpose of state estimation and control of freeway traffic. REFERENCES [1] C. Setchell and E. Dagless, “Vision-based road-traffic monitoring sensor,” IEE Proc.-Vis. Image Signal Process., vol. 148, no. 1, pp. 78–84, 2001. [2] Y. Wang, M. Papageorgiou, and A. Messmer, “A real-time freeway network traffic surveillance tool,” IEEE Trans. Control Syst. Technol., vol. 14, no. 1, pp. 18–32, Jan. 2006. [3] R. Bolla, F. Davoli, and M. Repetto, “On-line estimation of road traffic conditions by using mobile cellular networks,” presented at the 8th World Congr. Intell. Transport Syst., Sydney, Australia, 2001. [4] R. Narasimhan and D. Cox, “Estimation of mobile speed and average received power in wireless systems using best basis method,” IEEE Trans. Commun., vol. 49, no. 12, pp. 2172–2183, Dec. 2001. [5] V. Astarita, R. Bertini, S. d’Elia, and G. Guido, “Motorway traffic parameter estimation from mobile phone counts,” Eur. J. Oper. Res., vol. 175, no. 3, pp. 1435–1446, 2006. [6] M. Papageorgiou, “Applications of automatic control concepts to traffic flow modeling and control,” in Lecture Notes in Control and Information Sciences, A. Balakrishnan and M. Thoma, Eds. Berlin, Germany: Springer-Verlag, 1983, vol. 50. [7] M. Papageorgiou, J.-M. Blossevilee, and H. Hadj-Salem, “Macroscopic modelling of traffic flow on the Boulevard Périphérique in Paris,” Transportation Res. B, vol. 23B, pp. 29–47, 1989. [8] M. Papageorgiou, J.-M. Blossevilee, and H. Hadj-Salem, “Modeling and real-time control of traffic flow on the southern part of Boulevard Périphérique in Paris,” Transportation Res. A, vol. 24A, no. 5, pp. 345–359, 1990.
ALESSANDRI et al.: MODELING AND IDENTIFICATION OF NONLINEAR DYNAMICS FOR FREEWAY TRAFFIC
[9] J. Park and I. Sandberg, “Universal approximation using radial-basisfunction networks,” Neural Computation, vol. 3, no. 2, pp. 246–257, 1991. [10] M. Leshno, V. Ya, A. Pinkus, and S. Schocken, “Multilayer feedforward networks with a nonpolynomial activation function can approximate any function,” Neural Netw., vol. 6, no. 6, pp. 861–867, 1993. [11] F. Girosi, “Regularization theory, radial basis functions and networks,” in From Statistics to Neural Networks. Theory and Pattern Recognition Applications Subseries F, Computer and Systems Sciences. New York: Springer-Verlag, 1994, pp. 166–187. [12] V. Ku˚rková, “Approximation of functions by perceptron networks with bounded number of hidden units,” Neural Netw., vol. 8, no. 5, pp. 745–750, 1995. [13] J. Sjoberg, Q. Zhang, L. Ljung, A. Benveniste, B. Deylon, P. Glorennec, H. Hjalmarsson, and A. Juditsky, “Nonlinear black-box models in system identification: a unified overview,” Automatica, vol. 31, no. 12, pp. 1691–1724, 1995. [14] P. Gipps, “A behavioural car-following model for computer simulation,” Transportation Res. B, vol. 15, no. 2, pp. 105–111, 1981. [15] P. Gipps, “A model for the structure of lane changing decision,” Transportation Res. B, vol. 20, no. 5, pp. 403–414, 1986.
959
[16] H. Payne, “Models of freeway traffic and control,” in Simulation Council Proc., 1971, vol. 1, no. 1, pp. 51–61. [17] G. Whitham, Linear and Nonlinear Waves. New York: Wiley, 1974. [18] A. Sen and M. Srivastava, Regression Analysis, Theory, Methods, and Applications. New York: Springer-Verlag, 1990. [19] S. Chen, S. Billings, and P. Grant, “Nonlinear system identification using neural networks,” Int. J. Control, vol. 51, no. 6, pp. 1191–1214, 1990. [20] J. Suykens, B. D. Moor, and J. Vandewalle, “Nonlinear system identification using neural state space models applicable to robust control design,” Int. J. Control, vol. 62, no. 1, pp. 129–152, 1995. [21] R. Grino, G. Cembrano, and C. Torras, “Nonlinear system identification using additive dynamic neural networks—Two on-line approaches,” IEEE Trans. Circuits Syst. I, Fundam. Theory Appl., vol. 47, no. 2, pp. 150–165, Feb. 2000. [22] S. Billings and Q. Zhu, “Model validation tests for multivariable nonlinear models including neural networks,” Int. J. Control, vol. 62, no. 4, pp. 749–766, 1995. [23] T. Coleman, M. Branch, and A. Grace, Optimization Toolbox. Natick, MA: The MathWorks, Inc., 1999.
IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 17, NO. 4, JULY 2009
Modeling and Identification of Nonlinear Dynamics for Freeway Traffic by Using Information From a Mobile Cellular Network Angelo Alessandri, Senior Member, IEEE, Raffaele Bolla, Member, IEEE, Mauro Gaggero, and Matteo Repetto
Abstract—The high coverage of the territory by cellular networks and the widespread diffusion of mobile terminals aboard vehicles allow one to collect information on the traffic behavior. The problem of selecting a dynamic model to describe the freeway traffic by using the information available from a wireless cellular network is addressed by assuming the distribution of mobile terminals aboard vehicles to be uniform along the carriageway. Two different nonlinear parametrized models of freeway traffic are investigated: the first is an extension to a well-established macroscopic model, while the second is based on a black-box approach and consists in using a neural network to approximate the traffic dynamics. The parameters of such models are identified off line by a least-squares technique. Traffic measurements obtained from a cellular network are employed to identify and validate the proposed models, as shown by means of simulations. Index Terms—Freeway traffic model, identification, least squares, macroscopic model, mobile cellular network, neural networks.
I. INTRODUCTION N THIS brief, Dynamic models to describe the traffic on a freeway stretch are proposed that are based on the information about density, mean speed, and flow of vehicles that can be provided by a wireless cellular network, with a limited use of traditional measurement devices such as loop detectors, radars, and television cameras (see, e.g., [1]). Since such sensors are very expensive, a reduced number of them is placed along the roads, usually at long distances, thus providing only local information, possibly affected by inhomogeneities [2]. Here we shall exploit the traffic information that can be obtained from mobile terminals aboard vehicles, i.e., phones, handhelds, and personal computers in ongoing communications with a base station for a voice call or a data transfer (we shall denote such terminals as “active”). Starting from the ideas reported in [3], a section between two standard sensors is replaced with a cell of a wireless network. The cellular network technologies allow one to have a large number of small cells, 2–4 km in diameter, as compared with much larger distances between successive sensors. Moreover, the diffusion of wireless communications makes it very likely for vehicles to take aboard mobile terminals. The main
I
Manuscript received July 11, 2008; revised November 13, 2008. Manuscript received in final form January 22, 2009. First published May 19, 2009; current version published June 24, 2009. Recommended by Associate Editor C. C. Lee. A. Alessandri and M. Gaggero are with the Department of Production Engineering, Thermoenergetics, and Mathematical Models (DIPTEM), University of Genoa, I-16129 Genoa, Italy (e-mail: [email protected]; [email protected]). R. Bolla and M. Repetto are with the Department of Communications, Computer and System Sciences (DIST), University of Genoa, I-16145 Genoa, Italy (e-mail: [email protected]; [email protected]). Digital Object Identifier 10.1109/TCST.2009.2014242
drawbacks of standard acquisition devices (i.e., locality of traffic measurements and high costs) are overcome by acquiring information from mobile devices. Using the cellular network, one can directly collect the measurements of density of vehicles with active terminals aboard and flow of such vehicles moving from a cell to the next one. Moreover, recently developed techniques (see, e.g., [4]) allow one to determine the speed of active terminals, and hence of vehicles having them aboard. It is worth noting that, if required, one can integrate the information obtained from a cellular network with those available from traditional traffic sensors. The idea of collecting traffic measurements from mobile terminals aboard vehicles is reported also in [5], where the problem of estimating density and mean speed of vehicles is addressed by using a dynamic description of the traffic based on a partial differential equation that accounts for the conservation of the density of vehicles. In this respect, our contribution is that of developing two parametrized models that include also the speed dynamics and focusing on the identification of their parameters. Under the assumption on the uniform distribution of mobile terminals inside each cell of a wireless network, we propose two different approaches to model the freeway traffic. The first is an extension of a well-known macroscopic traffic model (see, e.g., [2], [6]–[8]). Such extension consists in adding a state variable that represents the percentage of vehicles with on-board active terminals with respect to the total number of vehicles in each cell. The second relies on the idea of using a neural network to approximate the traffic dynamics according to the black-box paradigm, and hence does not require all the fluid hypotheses used to construct the first. Such an approach is motivated by the excellent approximating properties of neural networks and the availability of efficient methods to tune the neural weights (see, e.g., [9]–[13]). Both the first and the second model are nonlinear in the state variables and depend nonlinearly on the corresponding parameters, which can be identified using least squares. In order to validate the models that result from the identification of the corresponding parameters, a simulation setup has been devised that is composed of: 1) a freeway traffic microscopic simulator that traces the behaviors of all the vehicles moving on the freeway (see [14], [15]) and 2) a cellular network simulator that takes into account voice calls or data transfers generated by on-board terminals. The rest of the brief is organized as follows. Section II describes the proposed extended macroscopic traffic model and the least-squares approach used for the identification of its parameters. In Section III, the same problems of modeling and identification are addressed via neural networks. Section IV reports the results obtained in a simulation case study. The conclusions are drawn in Section V.
1063-6536/$25.00 © 2009 IEEE
ALESSANDRI et al.: MODELING AND IDENTIFICATION OF NONLINEAR DYNAMICS FOR FREEWAY TRAFFIC
953
Fig. 1. Macroscopic model of traffic in a freeway trunk.
II. EXTENDED MACROSCOPIC TRAFFIC MODEL In the literature, there exist two main families of models to represent the traffic behavior along a road: the macroscopic ones, which deal only with aggregate quantities (e.g., density and mean speed of vehicles), and the microscopic ones, which describe the behavior of all the single vehicles. Using a microscopic model, one can easily simulate particular events such as slowing-downs, accidents, and carriageway restrictions. A microscopic simulator may be very detailed and precise, but requires a great computational effort. By contrast, in a macroscopic framework, the traffic is modeled like a fluid by using continuous variables, thus turning in less flexibility for simulating specific situations like the aforementioned ones; anyway, as compared with the microscopic approach, the computational burden is reduced. A macroscopic model of the traffic on a freeway trunk can be obtained by dividing the stretch into various sections, each described by state variables such as density and mean speed of vehicles (see Fig. 1). Moreover, one has to take into account the traffic flows leaving a section and entering the following one, as well as the input and output flows at the on- and offramps, respectively. We assume to perform a discretization at time instants , where and is the sampling period. At time , the following variables are defined in a freeway stretch made up of segments: 1) is the density of vehicles in section (in veh/km), ; 2) is the mean speed of vehicles in section (in km/h), ; 3) is the flow of vehicles leaving section and entering section (in veh/h), ; 4) and are the flows of vehicles entering and coming out of section , by means of on- and off-ramps (in veh/h), respectively; and are measured via standard sensors located at the ramps and are identically null if there exist no ramps in section . Second-order models are usually considered to describe the dynamics of freeway traffic by using such aggregate variables. Here we shall adopt the Payne–Whitham model [16], [17], by referring to the well-established approach developed by Papageorgiou and coworkers (see [2], [6]–[8], and the references therein). The density of vehicles in section varies from a time interval to the next, depending on the flows of vehicles that entered and left the section, according to the following balance equation:
(1)
, and is the length of section (in km), where which is assumed to be known. For the first and last sections, state-augmentation equations are considered, i.e., and , where the scalar noises and model the uncertainties in the densities at the boundaries. Concerning the speed dynamics, the following nonlinear equation is adopted:
(2) where and is a noise. The speed dynamics in the boundary segments is described by and , where and are scalar noises. Moreover, , and are positive constants and is the function representing the static speed-density relationship, where is the free speed, and is the critical density. Finally, the flow of vehicles between section and section is given by (3) and is a constant parameter. where The model given by (1), (2), and (3) can be extended by using the traffic measurements from a cellular network that refer to a subset of vehicles (i.e., those with an active terminal aboard). Such information is basically the density of active terminals in each cell of the network and the handoff flows. Moreover, methodologies are now available to estimate the speed of a moving active terminal by using various techniques (see [4] and the references therein). Thus, we shall present a model that allows one to predict the behavior of the overall traffic in a freeway stretch by using the information obtained from a cellular network. Such information refers to a subset of vehicles (i.e., density, mean speed, and flow only of the vehicles with an active terminal aboard). If just one mobile terminal were on each vehicle and one would know the total number of such terminals in each cell of the wireless network, the density of vehicles on each cell could be computed exactly. Unfortunately, the cellular networks do not know the precise cell inside which each terminal resides. As a matter of fact, in order to limit the signalling between the
954
IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 17, NO. 4, JULY 2009
Fig. 2. Extended macroscopic model of traffic in a freeway trunk.
mobile terminals and the network, the positions of the terminals are updated only with respect to the location area, which is a set of quite a large number of cells that represents the most precise information about the positions of the terminals. However, when a mobile terminal is active (i.e., when it is communicating with a base station), the cell on which it resides is known. In other words, a cellular network can supply information about the number of active terminals on each cell and the flows between adjacent cells due to handoffs that occur because of the movement of active terminals from one cell to the next one along the freeway. In the following, we denote the vehicles with an active terminal aboard as the “observed” ones. The freeway sections are replaced with network cells. The ratio between the number of vehicles with active terminals aboard and the total number of vehicles in cell at time is denoted by and regarded as a state variable for each cell together with density and mean speed. The resulting model is depicted in Fig. 2. To sum up, the description of the proposed model requires the definition of the following additional variables at time : 1) is the percentage of observed vehicles with respect to the total number of vehicles in cell ; 2) is the density of observed vehicles in cell (in veh/ km), ; 3) is the density variation of observed vehicles in cell (in veh/km) due to the entering/leaving from the ramps or the start/end of voice calls or data transfers, ; 4) is the flow of observed vehicles leaving cell and entering cell (in veh/h), . The balance of observed vehicles entering and leaving cell yields the following relationship:
where
, and (6)
, and is given by (3). Thus, where can be computed by means of (4) and (6) as follows:
(7) where . For the boundary segments, the percentages of observed vehicles are assumed to satisfy the following relationships: and , where and are scalar noises. To derive a dynamic equation for that accounts for the measurements provided by the cellular network, we substitute (6) into (1) so as to obtain
(8) where . We recall that the measures of the variables , and are available from the cellular network, as well as and can be easily determined by standard sensors located at the on- and off-ramps, respectively. Moreover, we assume to measure the mean speed of vehicles by averaging the speeds of active terminals; we denote these measures by . Likewise in [5], we suppose that the mean speeds of the observed vehicles and of all the vehicles approximately agree under the uniformity assumption, i.e., (9)
where dynamics for
. Using this equation and (1), we derive the , i.e.,
(4) where . If we assume the distribution of vehicles with active terminals to be uniform inside each cell, the following relationships hold: (5)
(note that in the following we shall introduce where measurement noises to model uncertainties in the measures). Therefore, keeping in mind (5), (6), and (9), let us define the measurement vector as
where
and, from now on, the notation stands for a generic column vector . In addition, let us introduce the following definitions at time : 1) is the vector of the state variables;
ALESSANDRI et al.: MODELING AND IDENTIFICATION OF NONLINEAR DYNAMICS FOR FREEWAY TRAFFIC
2) is the vector that collects all the uncertainties affecting the dynamics; 3) is the vector of the density variations of observed vehicles due to their entering/leaving from the ramps or the start/end of voice calls or data transfers; 4) and are the vectors containing the flows of vehicles at the on- and off-ramps, respectively. For the sake of compactness, let us denote the difference equation and the output mapping that represent the dynamic equations [i.e., (8), (2), and (7)] and measurement equations [i.e., (5), (6), and (9)] by means of the functions and , respectively: (10a) (10b) where is the vector of the parameters to be determined together with , and is the vector of the measurement noises. The model (10) depends nonlinearly on the parameters and . According to a classical least-squares approach (see, e.g., [18]), the noise vectors and are supposed to be “small,” so one can try to estimate and by minimizing a regression cost function based on a batch of data given by the measures , and for . On the basis of the foregoing, after choosing , the predicted state is given by
955
rameters (i.e., neural weights) to ensure a fixed approximation accuracy, especially in high-dimensional settings. More specifically, one-hidden-layer sigmoidal neural networks and radialbasis-function networks with tunable external and internal parameters enjoy the property of uniform approximation and guarantee a precision that grows at most polynomially with the dimension of the input of the function to be approximated (see, e.g., [9]–[12] and the references therein). It is worth noting that such a property does not hold for linear approximators. Another crucial aspect that may suggest the use of neural approximators is the possibility of using efficient methodologies for optimizing the neural weights (e.g., the Levenberg–Marquardt or Newton algorithms [13]). For the purpose of comparison with the extended macroscopic model proposed in Section II, we consider a model of the system (10) without all the fluid assumptions introduced to derive the function . Thus, an alternative expression of is searched for by using a neural mapping and the available measurements given by . Clearly, this entails the necessity of relying on the assumptions used to construct the function . Toward this end, according to well-established approaches available in the literature (see, among others, [19]–[21]), let us consider the black-box parametrized model
where
and the function
is a neural mapping that is parametrized by the -dimensional neural weights vector with the compact set . As in Section II, let us introduce a least-squares cost function that depends on and , i.e.,
(11) where and is some prediction of the initial state. A least-squares quadratic cost function that measures the quality of the fitting is
(12) where is a symmetric, positive definite matrix. An estimate of and can be obtained by minimizing the cost (12) under the constraints regarding the dynamic equations and the ranges of values that can be taken on by the parameters and the state variables. Thus, the problem reduces to search for and that minimize (12) under the constraints (11), . III. BLACK-BOX MODELING USING NEURAL NETWORKS The use of neural networks for the purpose of modeling is motivated by the fact that continuous functions can be approximated arbitrarily well on a compact set by large classes of neural networks with one or more hidden layer(s). Such a capability is referred to as “universal approximation property.” Besides, one-hidden-layer neural networks exhibit another powerful feature that consists in requiring a small number of pa-
(13) where is a symmetric, positive definite matrix. The optimal weights and parameter can be estimated by minimizing the cost (13) under the same type of constraints considered in the minimization of the cost (12) for the identification of the extended macroscopic model. IV. SIMULATION RESULTS A. Description of the Simulation Setup In order to evaluate the effectiveness of the proposed models, a simulation setup has been devised that is made up of a microscopic freeway traffic simulator and a cellular network simulator. The microscopic traffic simulator is based on the car-following model proposed by Gipps and adapted to the freeway environment [14], [15]. The Gipps model allows one to describe the typical driver behavior according to the physical laws of motion by taking into account the traffic conditions around each vehicle. The main idea behind such a modeling approach is that each driver reacts to a stimulus coming from other vehicles on the road, especially the previous and the following ones. Different situations of possible interest can be easily simulated, like, for example, accidents and lane restrictions, which are
956
IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 17, NO. 4, JULY 2009
Fig. 3. Estimates of the parameters of the extended macroscopic model during the optimization process.
Fig. 4. Behaviors of the state variables of the extended macroscopic model compared with the corresponding ones given by the microscopic traffic and cellular network simulators.
particularly well-suited to testing our models. We used the microscopic traffic simulator to reproduce the real behavior of vehicles on a freeway. The data generated by such simulator together with the cellular network simulator are used as measurements to compute the state variables of the macroscopic and neural models. The values of the percentages of observed vehicles, density, and mean speed of all vehicles given by such simulators are used for comparison with the same variables obtained by the proposed models.
As to the cellular network simulator, we considered the typical operations of a generic wireless network at the connection level since no details concerning the protocol (e.g., ETACS, GSM, UMTS, WiFi, WiMax) and the switching technique (circuit or packet) are taken into account. Radio resources are limited for each cell of the network, and hence the number of concurrent connections is limited as well. In the following, we shall focus on a simulation case study consisting of a one-way freeway stretch of 24 km covered by
ALESSANDRI et al.: MODELING AND IDENTIFICATION OF NONLINEAR DYNAMICS FOR FREEWAY TRAFFIC
Fig. 5. Autocorrelation functions of the residuals of q
957
(t); i = 4; 5; 6, obtained with the extended macroscopic model.
Fig. 6. Behaviors of some state variables of the neural black-box model compared with the corresponding ones given by the microscopic traffic and cellular network simulators.
12 cells of a cellular network (numbered sequentially from 0 to 11). Each cell is 2 km long, and there exist on- and off-ramps in cells 2, 6, and 10. The reference speed of each vehicle, the rate of entrance from the beginning of the road, as well as the rates of entrance and exit from the ramps, were randomly chosen according to uniform distributions in the ranges [70, 100] km/h, [500, 1000] veh/h, [230, 270] veh/h, and [170, 270] veh/h, respectively. The connections to the network were generated according to a Poisson distribution with mean proportional to the vehicle density [the proportional coefficient was chosen to be equal to 150 Erlang km/(s veh)]. The sampling period was taken equal to 10 s and each simulation was 2 h long.
For a quantitative evaluation of the model fitting, we applied standard autocorrelation tests on the residuals between the measured variables and their predictions given by the outputs of a model (see [22] for an introduction). Toward this end, let us define the residuals of a scalar measurement sequence , as the difference between and its corresponding prediction . The autocorrelation function of the residuals , is defined as
958
IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY, VOL. 17, NO. 4, JULY 2009
Fig. 7. Autocorrelation functions of the residuals of q
where
and .
(t); i = 4; 5; 6, obtained with the neural black-box model.
is the sample mean, i.e.,
B. Identification and Validation of the Extended Macroscopic Model In order to identify the parameters of the extended macroscopic model, we used the routine fmincon of the Matlab Optimization Toolbox (see [23] for details) to minimize the cost (12) with equal to and a batch of data that takes into consideration various events occurring on a freeway stretch (e.g., slowing-downs and accidents). The estimates of the parameters at each iteration of the optimization are shown in Fig. 3. The parameters converged asymptotically to the following values: 0.0031 h, 40.56 km /h, 47.97 veh/km, 105.14 km/h, 70.48 veh/km, and . For the sake of brevity, we report only the simulation results with the behaviors of the state variables of the identified model in the central cells (i.e., cells 4–6), with a slowing-down occurring in cell 4 (see Fig. 4). The dashed lines represent the “real-world values” of the state variables given by the microscopic traffic and cellular network simulators, while the continuous lines represent the state variables of the extended macroscopic model with the identified parameters. In Fig. 5, the autocorrelation functions of the residuals of , are shown, where the horizontal lines bound the 95% confidence bands. The autocorrelation functions lie outside the 95% confidence bands for a small number of lags, thus showing a satisfying behavior of the identified model. C. Identification and Validation of the Neural Black-Box Model We used one-hidden-layer feedforward neural networks with different numbers of sigmoidal units in the hidden layer to identify only the dynamics of the central cells (i.e., cells 4–6). Therefore, the neural networks had only 21 inputs and 9 outputs, which considerably reduced the computational burden for the selection of the weights and parameter . The minimization of the cost (13) with equal to was performed using the routine fmincon of the Matlab Optimization Toolbox. Among the various attempts, Fig. 6 shows the best results obtained by a one-hidden-layer feedforward neural network with 10 sigmoidal units in the hidden layer. The value of after the identification turned out to be equal to 0.8. By comparing such results with those obtained by the extended macroscopic
model (see Fig. 4), one can see that the speed dynamics is identified very well, while the percentage of on-board active terminals suffer from the problem of a small bias. As to the density of vehicles, the performances of the two models are quite similar. However, the autocorrelation functions of the residuals of , given by the neural black-box model lie outside the 95% confidence bands for a larger number of lags with respect to those of the extended macroscopic model, as shown in Fig. 7. V. CONCLUSION The problem of modeling and identifying the dynamics of freeway traffic by using the information available from active mobile terminals aboard vehicles has been faced according to two different paradigms. The first results from an extension to a well-known macroscopic model of freeway traffic. The second is based on the idea of using a black-box neural model, thus avoiding a lot of fluid assumptions that are required by the first approach. The analysis of the results obtained by the identification shows that the former approach has turned out to perform better than the latter on the overall. Future research will concern the use of the proposed models for the purpose of state estimation and control of freeway traffic. REFERENCES [1] C. Setchell and E. Dagless, “Vision-based road-traffic monitoring sensor,” IEE Proc.-Vis. Image Signal Process., vol. 148, no. 1, pp. 78–84, 2001. [2] Y. Wang, M. Papageorgiou, and A. Messmer, “A real-time freeway network traffic surveillance tool,” IEEE Trans. Control Syst. Technol., vol. 14, no. 1, pp. 18–32, Jan. 2006. [3] R. Bolla, F. Davoli, and M. Repetto, “On-line estimation of road traffic conditions by using mobile cellular networks,” presented at the 8th World Congr. Intell. Transport Syst., Sydney, Australia, 2001. [4] R. Narasimhan and D. Cox, “Estimation of mobile speed and average received power in wireless systems using best basis method,” IEEE Trans. Commun., vol. 49, no. 12, pp. 2172–2183, Dec. 2001. [5] V. Astarita, R. Bertini, S. d’Elia, and G. Guido, “Motorway traffic parameter estimation from mobile phone counts,” Eur. J. Oper. Res., vol. 175, no. 3, pp. 1435–1446, 2006. [6] M. Papageorgiou, “Applications of automatic control concepts to traffic flow modeling and control,” in Lecture Notes in Control and Information Sciences, A. Balakrishnan and M. Thoma, Eds. Berlin, Germany: Springer-Verlag, 1983, vol. 50. [7] M. Papageorgiou, J.-M. Blossevilee, and H. Hadj-Salem, “Macroscopic modelling of traffic flow on the Boulevard Périphérique in Paris,” Transportation Res. B, vol. 23B, pp. 29–47, 1989. [8] M. Papageorgiou, J.-M. Blossevilee, and H. Hadj-Salem, “Modeling and real-time control of traffic flow on the southern part of Boulevard Périphérique in Paris,” Transportation Res. A, vol. 24A, no. 5, pp. 345–359, 1990.
ALESSANDRI et al.: MODELING AND IDENTIFICATION OF NONLINEAR DYNAMICS FOR FREEWAY TRAFFIC
[9] J. Park and I. Sandberg, “Universal approximation using radial-basisfunction networks,” Neural Computation, vol. 3, no. 2, pp. 246–257, 1991. [10] M. Leshno, V. Ya, A. Pinkus, and S. Schocken, “Multilayer feedforward networks with a nonpolynomial activation function can approximate any function,” Neural Netw., vol. 6, no. 6, pp. 861–867, 1993. [11] F. Girosi, “Regularization theory, radial basis functions and networks,” in From Statistics to Neural Networks. Theory and Pattern Recognition Applications Subseries F, Computer and Systems Sciences. New York: Springer-Verlag, 1994, pp. 166–187. [12] V. Ku˚rková, “Approximation of functions by perceptron networks with bounded number of hidden units,” Neural Netw., vol. 8, no. 5, pp. 745–750, 1995. [13] J. Sjoberg, Q. Zhang, L. Ljung, A. Benveniste, B. Deylon, P. Glorennec, H. Hjalmarsson, and A. Juditsky, “Nonlinear black-box models in system identification: a unified overview,” Automatica, vol. 31, no. 12, pp. 1691–1724, 1995. [14] P. Gipps, “A behavioural car-following model for computer simulation,” Transportation Res. B, vol. 15, no. 2, pp. 105–111, 1981. [15] P. Gipps, “A model for the structure of lane changing decision,” Transportation Res. B, vol. 20, no. 5, pp. 403–414, 1986.
959
[16] H. Payne, “Models of freeway traffic and control,” in Simulation Council Proc., 1971, vol. 1, no. 1, pp. 51–61. [17] G. Whitham, Linear and Nonlinear Waves. New York: Wiley, 1974. [18] A. Sen and M. Srivastava, Regression Analysis, Theory, Methods, and Applications. New York: Springer-Verlag, 1990. [19] S. Chen, S. Billings, and P. Grant, “Nonlinear system identification using neural networks,” Int. J. Control, vol. 51, no. 6, pp. 1191–1214, 1990. [20] J. Suykens, B. D. Moor, and J. Vandewalle, “Nonlinear system identification using neural state space models applicable to robust control design,” Int. J. Control, vol. 62, no. 1, pp. 129–152, 1995. [21] R. Grino, G. Cembrano, and C. Torras, “Nonlinear system identification using additive dynamic neural networks—Two on-line approaches,” IEEE Trans. Circuits Syst. I, Fundam. Theory Appl., vol. 47, no. 2, pp. 150–165, Feb. 2000. [22] S. Billings and Q. Zhu, “Model validation tests for multivariable nonlinear models including neural networks,” Int. J. Control, vol. 62, no. 4, pp. 749–766, 1995. [23] T. Coleman, M. Branch, and A. Grace, Optimization Toolbox. Natick, MA: The MathWorks, Inc., 1999.
