Manifold Learning Method with OD Matrix as Latent
Variables and Surrogate Modeling for Pedestrian
Network Design
INU021002

Mizuki Ogawa
ogawa@bin.t.u-tokyo.ac.jp
Eiji Hato
hato@bin.t.u-tokyo.ac.jp
Kenta Ishii
kenta.ishii.py@hitachi.com

Abstract
This paper focuses on pedestrian behavior and attempts to construct a method to estimate OriginDestination traﬃc distribution and route choice parameters simultaneously. In the previous estimation
methods, it is not easy to guarantee the uniqueness of the estimated solution because the objective function
is set individually. In this paper, we propose an estimation method that guarantees the uniqueness ofthe solution by ﬁnding the nearest point between the observed manifold and the model manifold byintroducing manifold learning, which is a ﬁeld of information geometry. Moreover, we conducted numericalexperiments to conﬁrm the stability of the solution using the proposed method and the previous estimation
method. In addition, we developed a surrogate model to speed up the Pareto frontier search, a set of
Pareto solutions, to deal with the sizeable computational load of the network design problem based on theestimated parameters. In a study of a park development plan for an actual site, we achieved a signiﬁcant
speedup over the usual computation time by having a feed-forward neural network learn network patterns
as input and objective function values as output.

Ogawa and Hato, and Ishii

1

INTRODUCTION
The prolonged period of COVID-19 is changing the way urban space is being used, as people arechanging their behaviors by restricting self-restraint movement(Sakai et al. (1)). Thus, while thereallocation of open spaces for distance is becoming more critical, there is an increasing need for
ﬂow prediction based on real-time measurements of human ﬂow and density.
In order to understand the actual state of movement in urban space and to clarify the effects of policies, it is necessary to describe dynamic activity patterns based on partial observation
information from various sensors(Hato (2),Oyama and Hato (3)).
At the same time, it is essential to estimate the OD matrix, which is an input variable forthe allocation model. Based on the accumulation of many studies on OD estimation, we propose a
new methodology for the optimal design of urban spaces by constructing a simultaneous estimation
method for the parameters of the behavioral model.
However, the parameters of the behavioral model and the OD matrix of pedestrians arenested in each other, and it is theoretically and allegorically challenging to estimate them uniquely.
In this study, we propose an estimation method that theoretically guarantees the uniqueness of thesolution by treating the OD matrix as a latent variable of the link ﬂow using the observation manifold and treating the behavioral model as a recursive model manifold. By learning the manifold,
the OD matrix and the sequential choice probabilities of links or cells are estimated simultaneously.
On the other hand, it is necessary to perform a multidimensional evaluation for a vast number of combinations to design cities to meet new needs by arranging open spaces and developing
networks. We can optimize through Pareto frontier analysis based on trade-off relationships among
policy variables. However, it causes a signiﬁcant computational load, and we need to speed up thecomputation speed. In recent years, there has been much interest in applying surrogate models inmolecular biology and the optimal design of materials. Recently, there has been much interest in
applying surrogate models in molecular biology and the optimal design of materials. In this paper,
we propose a surrogate model, a new method to achieve faster computation speed by replacing
the input and output of a theoretical model with an approximate model through machine learning.
In this study, we aim to achieve fast policy simulation by transferring the pedestrian allocation
calculation results to a neural network model to calculate approximate solutions for many policy
simulation results.
In this study, the Chapter2 reviews the simultaneous estimation method of OD matrix androute choice models and the method of models with latent variables. We also review the applicationof surrogate models in engineering. In Chapter3, we explain the research framework and network
notation, and manifolds. In Chapter4, we construct a method for simultaneous estimation of OD
matrix and route choice models using manifold learning and discuss the properties of the solution.
Finally, in Chapter5, we conducted numerical experiments to conﬁrm the method’s performance,
followed by the estimation using real data in Chapter6. Chapter 7 explains the network design
method using a surrogate model and shows its speedup. Finally, in Chapter 8, conclusions andfuture research.
LITERATURE REVIEWIn this chapter, we summarize the studies that have been proposed on methods for simultaneous
estimation of OD matrix and route choice models and estimation methods for models with latentvariables. We also review surrogate models used as alternatives to simulation and consider using
them in this research.

Ogawa and Hato, and Ishii

2

Simultaneous estimation of OD matrix and route choice models
Liu and Fricker (4)constructed a method for simultaneous estimation from link ﬂow only by dividing OD matrix estimation and link dispersion parameter estimation into two stages based on
Dial allocation(Dial and Robert (5)).In this study, an evaluation index of the estimated values is
the mean square error between the observed link trafﬁc and the link trafﬁc calculated using theestimated values. However, this index cannot evaluate the estimation error across the OD matrix
and link cost parameters because it does not evaluate the trafﬁc dimension of each OD. Yang et al.
(6), Balakrishna et al. (7), and Wang et al. (8) proposed a simultaneous estimation method using
generalized least squares (GLS) as the objective function, underpinned by stochastic equilibrium
allocation for Dial NW. Moreover, Lo and Chan (9) proposed a simultaneous estimation method
using Generalized Least Squares (GLS) as the objective function, assuming that OD matrix follows a normal distribution. However, these estimation methods need to prepare the information
observed not only link ﬂow but also OD matrix.
Based on these studies, it is necessary to develop an estimation method that can guarantee
the uniqueness of the solution in the NW where its OD matrix is unknown. In addition, no research
has proposed a simultaneous estimation method within an acceptable framework for cyclic routes
to represent pedestrian stay and migration. Therefore, it is necessary to allow Markov chain allocation within the framework of simultaneous estimation methods for the analysis of pedestrian behavior. This study attempts to construct a simultaneous estimation method using the Oyama and Hato(10) that can express cyclic paths. In the route choice model, the Oyama and Hato (10) can represent the cyclic path by inserting a time-structure NW on the Recursive Logit model proposed by
Fosgerau et al. (11), which assumes Markov chain allocation. This study attempts to construct asimultaneous estimation method using the Oyama and Hato (10) that can express cyclic paths.
Model estimation with latent variables
We can obtain the link ﬂow using the current observation technology. However, we cannot observe
the origin and destination points of the observed link ﬂow due to privacy concerns. Therefore,
when prior information on OD matrix is not available, simultaneous estimation of OD matrix androute choice model can be regarded as a mixture distribution estimation problem with OD matrix
distribution as a latent variable.
Then, in this section, we review the estimation problem of the model with latent variables.
In the transportation ﬁeld, Bhat (12) and Vij et al. (13) estimated models for travel demand estimation and transportation choice, respectively, using the Expectation and Maximization (EM)algorithm to handle individual and group heterogeneity. In the observation error problem in route
choice models, Oyama and Hato (3) dealt with the estimation problem using the links of the observed data as latent variables. We can interpret this study as a method that assumes that in the E
step of the expected value calculation in the EM algorithm, we always belong to the class with themost signiﬁcant probability and then perform likelihood maximization in the subsequent M step.
Thus, the EM algorithm is used for the model estimation with latent variables, but one of
the challenges of this method is local solutions in the solution. For this problem, an approach using
a manifold concept has been proposed from the viewpoint of information geometry Akaho (14).
The set S = { f (x|θ )}of random variables with parameter θ can be viewed as a space (manifold) with θ as the local coordinate system. This space does not have a distance structure like Euclidean space, and space may be curved. However, it can be regarded as locally Euclidean space,
and the local coordinates θ reduce the dimension and represent the curved space. In this way, byOgawa and Hato, and Ishii

3

FIGURE 1 Decomposition of observed link trafﬁc by the latent variableﬁtting a local coordinate system to a high-dimensional space and reducing the number of dimensions, we can ﬁnd a model manifold that best ﬁts the observed data called manifold learning. If weconsider S as a space of probability distributions, called an exponential family of distributions, adual coordinate system, the e-coordinate system θ and the m-coordinate system η , is introduced.
The observed data can be represented as a single point in the η -coordinate, and the projection of
the dual coordinate system θ onto the model space M is orthogonal to the θ -coordinate. This iscalled the m-projection, and the projection from the model space to the observation distribution iscalled the e-projection. There is a deep relationship between divergence, representing the gap between probability distributions and e-projection and m-projection. It is known that e-projection andm-projection between sub-spaces of an exponential distribution family give the minimum value of
divergence.
Based on the above characteristics of e-projection and m-projection, Csiszaér and Tusnaédy(15) deﬁned the model manifold as the joint distribution model p(x, z) of the latent variable z andthe observed variable x with the model parameter ξ as the coordinate system, and the observation
manifold as the product q̂(x)q(z|x) of the empirical distribution q̂(x) of the observed variable obtained as observation and the parameter q(z|x).They expressed the degree of separation between
each manifold by KL divergence, which indicated the similarity between probability distributions.
They proposed the em algorithm, which aims at minimizing KL divergence by alternately projecting between each manifold(Figure ??). Moreover, Amari (16) showed that the em algorithm isequivalent to the EM algorithm when dealing with probability distribution spaces and that the em
algorithm has a uniqueness of solution if there exists a duality relation where the observed data
is a mixture distribution family. The estimated model is an exponential distribution family. In
the transportation ﬁeld, Fosgerau et al. (17) uses the duality of log-sum variables to generalize themodel to reﬂect exploratory behavioral norms.

Ogawa and Hato, and Ishii

4

FIGURE 2 e-projection and m-projection between model manifold and observed manifold

From the above, we can understand the link ﬂow of pedestrian by the current technology.
However, it is difﬁcult to understand the OD matrix, and there are few studies to predict the pedestrian scale OD matrix. Therefore, we developed a method for simultaneous estimation of OD
matrix and route choice model using only observed link ﬂow for a route choice model representing
pedestrian. In formulating OD matrix as a latent variable, we propose a method that can guarantee
the uniqueness of the solution based on manifold learning of information geometry.
Speeding up by surrogate models
The surrogate model can simulate physical phenomena as a surrogate model, and it can replacecomputationally extensive simulations such as stress analysis. In stress analysis, we can understandthe performance of structure by inputting the geometry, decomposing the mesh, and analyzing thestress in mesh units. When we increase the mesh resolution, we can analyze with high accuracy,
but the computational scale becomes very large. Thelin et al. (18) has developed a prediction
model that uses a surrogate model to replace the FEA portion of the fatigue life simulation, which
is necessary for determining jet engine components. A model that can be analyzed in real-timehas been constructed by learning ﬁve shape parameters as input and the results of FEA analysis
as output. Mai et al. (19) and Chang and Cheng (20) also showed that a neural network surrogate
model could replace the analysis of nonlinear structures using ﬁnite element analysis, signiﬁcantly
reducing computational cost and guaranteeing convergence.
Surrogate models are adequate for computationally demanding and complex simulations.
The network design problem(Hao (21),Farahani et al. (22),Oyama and Hato (23)) we attempt is
to ﬁnd an optimal multi-purpose solution with a large number of network patterns. The allocation
computation of the Recursive Logit model(Fosgerau et al. (11))used in this research has a recursive
computation process, which causes a huge computational load. Therefore, we propose a surrogate
model using a feed-forward neural network to replace the allocation calculation.
FRAMEWORK AND NOTATIONSIn this chapter, we explain about research framework and deﬁne the variables to be used in themodel, and from the viewpoint of information geometry, we deﬁne the model manifold and theobservation manifold.

Ogawa and Hato, and Ishii

5

Research framework
We introduce the framework of this research(3). In this study, we construct a method for simultaneous estimation of OD trafﬁc and route choice models and design a network using the obtained
parameters. Then, we constructed a surrogate model to reduce the computational load of calculating enormous network patterns and allocations and accelerate the computation.

Estimation of Parameter and OD matrix
OD matrix

Dynamic Route Choicelatent valuables
O

D

Parameter
OD pair

OD matrix and Route choice paremeterby Manifold learning

Network Design
Train data

Surrogate modeling
Input
Network pattern
OutputObjective function values
ex. expected utility, trip length

Multi-purpose
optimal solution
ParateFrontier

Urban design senario
Speed up allocation calculations.

FIGURE 3 Research Framework : Simultaneous estimation of OD matrix and route choiceparameters and Network Design using Surrogate model
Network notation
We deﬁne a network as a directed graph G = (A, E) where the node-set is A and the link set is E
(Figure 4).We deﬁne a network as a directed graph G = (A, E) where the node-set is A and the link
set is E.For a connected node k ∈ A, we deﬁne A(k) to be the subset of nodes connected to k, andl ∈ E to be the link connecting nodes k and a.And we deﬁne the deﬁnite term of the transition utility

Ogawa and Hato, and Ishii

6

from k to a ∈ A(k) as v(a|k), and its explanatory variable vector as xl . In this paper, we describe
the network without time properties for the sake of generalization, but It can also be applied totime-structured networks such as in Oyama and Hato (10).

FIGURE 4 Network representation:link l connects node k and node a ∈ A(k)Model Manifold
We formulate the link ﬂow derived from each OD using a route choice model. Our goal is toconstruct a simultaneous estimation of OD matrix and route choice model that can also represent
the migratory behavior of pedestrians. Therefore, we use the RL model (Fosgerau et al. (11)),which calculates the link choice probability based on the Markov chain allocation that can express
cyclic ﬂows, to formulate the allocated trafﬁc volume on a logit base.
Let z be a set of OD pairs. We assume that a traveler moves between some OD pairsz j ∈ z by choosing nodes sequentially from the origin node o(z j ) to the destination node d(z j ) in a
directed graph G = (A, E). Let u(a|k) = v(a|k) + µε (a) be the instantaneous utility that the traveler
can obtain when moving from the current node k to the neighboring node a ∈ A(k). The error term
ε (a) is assumed to be an i.i.d extreme value distribution of type I with zero mean, and µ is a scaleparameter. In this case, the transition probability from node k to a is given by the MNL model asfollows(1).
exp(v(a|k) +V d (z j )(a))
p(a|k, z j , θ ) =(1)∑a′ ∈A exp(v(a′ |k) +V d (z j )(a′ )θ is the vector of explanatory variable parameters, and v(a|k) = θ T xl . xl is explanatory variable
vector of link l. The state value function V d (z j )(k), which is the expected downstream utility, canbe formulated recursively using the Bellman equation as follows(2).
V d (z j )(k) = E[max v(a|k) +V d (z j )(a) + µε (a)](2)
aAnd from the assumption of the error term, equation (2) can be transformed into a log-sum form
as follows (3) .
{d(z j )
1
µ ln ∑a′ ∈A(k) e µ {v(a|k)+V (a)} (k ̸= d(z j ))d(z j )
V
(k) =
(3)
0
(k = d(z j ))

Ogawa and Hato, and Ishii

7

Using the RL model, we calculate the trafﬁc volume vl j of link l derived from OD pair z j . Let
σ (l, j) denote the set of paths from o(z j ) to the starting node k of link l. Each path is deﬁned asa set of node sequences {sr0 , sr1 , ..., srI , sri ∈ A(sri−1 )} = r ∈ σ (l, j) where srI = k. Let q j be the OD
matrix of OD pair z j . Then vl, j can be expressed as follows(4).
vl, j =q j p(a|k, z j ; θ )I−1∑ ∏ p(sri+1|sri , z j ; θ )r∈σ (l, j) i=0

I−1
1=q j exp( {v(a|k) +V d(z j ) (a) −V d(z j ) (o) + ln ∑ ∏ p(sri+1 |sri , z j ; θ )})
(4)
µ
i=0
r∈σ (l, j)
The fourth term on the right-hand side of (4) is a log-sum value that represents the expected valueof the maximum utility obtained from o(z j ) to node k, and can be transformed as follows(5).

1vl, j = q j exp( {v(a|k) +V d(z j ) (a) + V̄ o(z j ) (k) −V d(z j ) (o)}(5)µ
V̄ o(z j ) (k) can be formulated recursively by applying the idea of (5) as follows(6).
{o(z j ) ′1′
µ ln ∑k′ ∈Ā(a) e µ {v(a|k )+V̄ (k )} (k ̸= o(z j )))o(z j )
V̄
(a) =
(6)
0
(k = o(z j ))Here, Ā(a) is deﬁned as the set of nodes adjacent to node a on the graph Ḡ that inverts the directed
graph G. We can express the link ﬂow for each OD using only structural and explanatory parameters. We can apply to Dial NW, which is consistent with Vliet (24) by this formulation. When weapply this formulation to NW with cyclic paths, we need a negative constraint in the determinantterm of utility. From the above, the simultaneous distribution of xl is explanatory variable vectorof link l xl and OD pair z j can be formulated as follows(7), using the parameter ξ = {θ , w}.
vl, j∑l ∑ j vl, j
wj
1
=
(7)
exp( {v(a|k) +V d(z j ) (a) + V̄ o(z j ) (k) −V d(z j ) (o)}B(ξ )µw is a weight parameter satisfying ∑ j w j = 1, and B(ξ ) is a normalization term, xl is explanatory
variable vector of link l and OD pair zl . (7) can be viewed as a manifold with ξ as its coordinate
system on space S = p(x, z) of the entire probability distribution, which we deﬁned as the model
manifold.
p(xl , z j ; ξ ) =Observed Manifold
In this study, we treat the link trafﬁc of travelers in the network as an observation. Since the linkﬂows vl can be assumed to follow the probability of occurrence q(x) of the explanatory variablesof the link, the observed link ﬂow can be expressed as q̂(xl ) = vBl . The B represents the sum of theobservable link trafﬁc on the network. In a static network, link ﬂow is the number of transitions
between points, but in a time-structured network, we can also consider the amount of stay at apoint as an observation.
When we use link trafﬁc as an observation, the problem is that we cannot observe the
OD of each traveler. Therefore, the empirical distribution obtained from the observation does
not correspond to the distribution model, which is the simultaneous distribution of OD matrix

Ogawa and Hato, and Ishii

8and link ﬂow, and the OD matrix distribution becomes a latent variable. Moreover, we cannot
uniquely deﬁne the observation points on space S = p(x, z) of the entire probability distribution. x
is explanatory variable vector and z is OD pair set. From the above, we deﬁne the manifold D byobservation as follows(8).
D = {q̂(x)q(z|x)}(8)q(z|x) is an arbitrary conditional distribution, deﬁned as a parameter corresponding to the coordinate systems z and x. Equation(8) means equivalent to deﬁning the observation as a mixture
distribution.
EM ALGORITHM AND STRUCTURE ESTIMATIONIn this chapter, we focus on the divergence between the model manifold and the observed manifold and the ﬂatness of the manifold and propose an estimation algorithm that combines structure
estimation and em algorithm to guarantee the uniqueness of the converged solution.
Objective function KL divergence
In the simultaneous estimation of the OD matrix and route choice model, we deﬁne the KL divergence as the objective function, representing the degree of separation between the two manifolds.
The KL divergence of two probability distributions f (x), g(x) is expressed as follows(9).
KL( f ||g) =∫f (x)[log f (x) − log g(x)]dx(9)

The KL divergence is non-negative and is not strictly a distance, but it is used as the degree ofseparation between probability distributions. When we do not know the actual probability distribution f (x), we need to estimate the actual probability distribution f (x) from the observed data.
In this case, we can estimate the model distribution by bringing it closer to the actual distribution,and the goal of estimation is to minimize the KL divergence. Based on the manifold deﬁned in theprevious section, we formulate the KL divergence(10).
q̂(x)q(z|x)
p(x, z|ξ )x zParameter ξ = {θ , w} is estimated by minimizing KL divergence as follows(11,12,13).
KL(q(z|x)||ξ ) = ∑ ∑ q̂(x)q(z|x) ln(10)θ̃ =θ KL(q(z|x)||ξ )w̃ =w KL(q(z|x)||ξ )
B
q̃ =
w̃
B(ξ̃ )

(11)
(12)
(13)

em algorithm
The em algorithm (Amari (16)) searches for the nearest point between the observation manifold
and the model manifold by alternating e-step and m-step. The e-step refers to the projection froma point on the model manifold onto the observation manifold (e-projection), and the m-step means
the projection from the observation manifold onto the model manifold (m-projection). In thissection, we will check the uniqueness of the solution at each step.

Ogawa and Hato, and Ishii

9e stepIn the e-step, we ﬁx ξ and consider the e-projection onto the observation manifold that minimizes
KL divergence with respect to q(z|x) Since q(z|x) is a probability distribution, it can be expressed as∑z j ∈z q(z j |x) = 1. In this case, we can express the Lagrangian of the KL divergence minimization
problem by using the undecided multiplier λl as follows(14).
q̂(xl )q(z j |xl )− ∑ λl ( ∑ q(z j |xl ) − 1)p(xl , z j |ξ )xl ∈x z j ∈zxl ∈xz j ∈zFrom KKT conditions, we can transform as follows(15).

L(q(z|x), λ |ξ ) = ∑ ∑ q̂(xl )q(z j |xl ) ln

(14)q(z j |xl )∂ L(q(z|x), λ |ξ )= q̂(xl ) ln
− ηl
(15)
∂ q(z j |xl )p(xl , z j |ξ )where ηl = q̂(xl ) ln q̂(xl ) + 1 − λl . We can transform equation (15) from Bayes’ theorem using therelation p(xl , z j |ξ ) = p(z j |xl , ξ )p(xl |ξ ) as follows(16).
q(z j |xl )∂ L(q(z|x), λ |ξ )(16)
= q̂(xl ) ln−τ
∂ q(z j |xl )p(z j |xl , ξ )where τl is a constant. In other words, when ξ is ﬁxed, the q(z j |xl ) that minimizes the KL divergence can be obtained as follows(17,18).
q(z j |xl ) =p(z j |xl , ξ )(17)p(xl , z j |ξ )=(18)
∑z j ∈z p(xl , z j |ξ )We can see that the e-projection is uniquely determined for the observational manifolds that wedeal with in this study.
m step
Next, ﬁx q(z|x) and consider the m-projection to the model manifold that minimizes the KL divergence with respect to ξ . When q(z|x) is ﬁxed, the ﬁrst term of the KL divergence ((10)) is constant,and minimization of KL divergence requires maximization of the second term, ln p(x, z|ξ ). Therefore, it is equivalent to considering the following maximization problem(19).
L(ξ |q(z|x)) = ∑ ∑ vl q(z j |xl ) ln p(xl , z j |ξ )(19)l∈E z j ∈zHere, the model manifold(7) contains structural parameters and is not an exponential family of
distributions. Therefore, there is no guarantee that the m-projection is uniquely determined, andthe uniqueness of the solution of the em algorithm cannot be guaranteed (Amari (16)).
4.3 Nested pseudo em algorithm
To guarantee that the em algorithm has a single solution, we propose the nested pseudo-em algorithm (NPem), which combines the nested pseudo-likelihood algorithm (NPL) (Aguirregabiria and Mira(25))with the em algorithm. The parameter ξ estimated by the em algorithm satisﬁes Equation(17),which can be transformed by substituting it into Equation(10) as follows(20).

Ogawa and Hato, and Ishii

10KL(q(z|x)||ξ ) = ∑ ∑ q̂(x)p(z|x, ξ ) lnxzq̂(x)p(z|x, ξ )p(x, z|ξ )= ∑ q̂(xl ) ln q̂(xl ) − ∑ q̂(xl ) ln p(xl |ξ )
l∈E

(20)
(21)

l∈E

As a result, since the structure parameters inside the em algorithm are ﬁxed in the Inner Algorithm,
the model manifold can be treated as an exponential distribution family, and the em algorithm
has a unique convergence point. If the solution converges, the convergence point is uniquelydetermined for the entire estimation algorithm. We will show the estimation ﬂow of the NPemalgorithm(Figure 5).

Ogawa and Hato, and Ishii

11

FIGURE 5 NPem algorithm is composed of Inner algorithm and Outer algorithm.

NUMERICAL EXPERIMENTIn this chapter, we conﬁrm the stability of the estimates of the proposed method by comparing with
Yang et al. (6) method of generalized least square (GLS) on the dial NW (Figure 6). In this study,
we do not assume equilibrium on the NW. In the experiment, we deﬁne the instantaneous utility

Ogawa and Hato, and Ishii

12obtained by a traveler’s node transition as follows(22).
v(a|k) = θlength LCl + θshop DCl(22)
where LCl is the link length of link l and DCl is the store dummy of the endpoint node a of linkl. The parameter θ = (−1.0, 2.0), and we set OD pair z1 = (1, 16), z2 = (5, 16), z3 = (2, 16), z4 =(1, 12) and OD matrix q1 = 200, q2 = 300, q3 = 300, q4 = 200.
Station(2, 0)

1(1, 0)(2, 0)

5(1, 0)(2, 0)

9(1, 0)

13(2, 0)(2, 0)

2(1, 0)(2, 1)

6(3, 0)

10(2, 0)(1, 0)

14

OriginNode(2, 0)(2, 0)

3

4(1, 1)(2, 0)

7Shop(1, 0)

11(2, 0)(3, 0)

15(2, 0)(1, 0)

8(1, 0)

12

Touristfacilities (1, 0)

16

Destination
Node

FIGURE 6 Dial NW for simulation. 3 origin and 2 destination.

We compare the mean, unbiased variance, 95% conﬁdence interval, minimum and maximum values, and RMSE of the parameter estimates for the estimation accuracy. The RMSE is used
as a parameter accurate value reproduction rate index, and the deﬁnition is as follows(23). Γ̂ is theactual parameter value, Γi is the ith estimated parameter, and n is the number of parameters.
√∑ni=1 (θ̂i − θi )2
RMSE(θ̂ ) =(23)n
We can see in Table 1 that there is no observation error in the link ﬂow and that the initial valuesof the parameters are generated with a uniform distribution for 100 trials. In the GLS, the mean of
θlength is −1.42, and there are cases where the minimum value is −8.9826, which is a signiﬁcantdecrease in the estimation accuracy. On the other hand, the average of θlength of the proposed
method is 1.0004, which means that the actual value can be estimated. RMSE, the proposed
method is 0.0910 compared to 1.0608 for GLS, indicating a high estimation accuracy.
Table 2 shows the estimation results when there is an observation error in the link trafﬁc.
Here, we assume that the observation error follows a normal distribution with variance σl = 5 forall links l. When observation errors exist, the RMSE of GLS is 8.0594, while that of NPem is6.3608, conﬁrming the superiority of the estimation accuracy.

Ogawa and Hato, and Ishii

13

TABLE 1 Estimation with simulated data in dial NW
Variable True value Method MeanSD
CI 95% Min
Max
Length-1.00
NP_em -1.0004 0.0018 0.0003 -1.0015 -0.9887
GLS-1.4200 1.5110 0.2961 -8.9826 -0.9994
Shop
2.00
NP_em 1.9980 0.0058 0.0011 1.9931 2.0200
GLS
1.9809 0.0644 0.0126 1.7634 2.0013
OD1200
NP_em 199.90 0.5076 0.0994 199.57 203.23
GLS
198.81 4.0614 0.7960 184.97 200.24
OD2300
NP_em 299.89 0.2494 0.0488 299.52 300.59
GLS
300.99 3.3748 0.6614 299.87 312.42
OD3300
NP_em 300.14 0.0452 0.0885 298.14 300.51
GLS
298.95 3.5661 0.6989 286.70 300.17
OD4200
NP_em 200.10 0.4229 0.0828 197.59 200.38
GLS
201.75 5.9140 1.1591 199.71 221.77
RMSE
NP_em 0.0910
GLS
1.0608TABLE 2 Estimation with simulated data in dial NW with observation error
Variable True value Method
Length
-1.00
NP_em
GLS
Shop
2.00
NP_em
GLS
OD1
200
NP_em
GLS
OD2
300
NP_em
GLS
OD3
300
NP_em
GLS
OD4
200
NP_em
GLS
RMSE
NP_em
GLS

Mean
-0.9729-1.3943
1.8713
1.9442203.77211.00286.05303.82306.05284.08199.35199.236.36088.0594SD
0.00281.46760.00580.05500.817718.2249
0.23442.70440.535716.79210.61697.2418

CI 95% Min
0.0005 -0.9752
0.2876 -12.4856
0.00111.86330.01081.75090.1603202.86
3.5720167.300.0459285.550.5301302.920.1050302.703.2912263.920.1209194.86
1.4194
194.77Max-0.9524-0.9908
1.90362.0015209.46222.93286.57316.18306.64312.25199.78228.89REAL DATA ESTIMATION
We estimate the actual data based on the proposed method. We use probe person data and mobile
spatial statistics data in the Toyosu area in Tokyo. In order to reduce the computational load, we
introduce a multi-scale network. In addition, we use a simultaneous estimation method that uses
link transition probability and link ﬂow as observations, and describe the details of the estimation
method in the appendix.
Data details and network settings
we will show the details of the data used in Tables 3.
We introduce a multiscale network consisting of grids of different sizes. As NW data, weconsider a multiscale network consisting of grids of different sizes. A node represents each grid,and links represent connections to neighboring grids. Each grid is connected if it is tangent to theOgawa and Hato, and Ishii

14

other grids, and the same applies to grids of different sizes (Figure 7). In this study, the explanatory
variables are the distance between cells, Store area, park area, number of high-rise buildings, anddummy variables for station existence, and the length of the road in the cell(Table 4).
TABLE 3 Prove person data and Mobile Spatial data
Prove person data
Survey Term
Survey Method
Survey Area
Number of Monitors
Number of Trips
Number of Valid Trips
Mobile spatial data
Survey Term
Survey Area
Monitors
Grid size
Time Unit

data detail30th Sep.2020 29th Nov.2020
Smart Phone Sensor
People living or working in Toyosu Area
296
37920
37439
data detail
1st Jul.2019 31th Jul.2020
All over the country
People use telephone of NTT docomo
500m mesh
1 hour(24Class)

TABLE 4 Explanatory variables
Explanatory Variables Overview
Distance between grids Shortest path length between the center of PP data in the adjacent grid
Store
Area of store in each 50m grid
Park
Area of park in each 50m gridHigh building
Number of building over 20th ﬂoor in each 50m grid
StationNumber of Station in each 50m gridRoad length
Road length in the cell divided by the cell area

Ogawa and Hato, and Ishii

15

2.5km

500m

100m

500m mesh

100m mesh

50m mesh

Connection between deﬀerent grid size

FIGURE 7 Multi-scale network and connection between different grid size

Estimation Result
We show in Table 5 the estimation results using the mobile spatial statistics data for each period andthe transition probabilities obtained from the PP data. The distance between the cells parameter isnegative in all periods, and the t-value is also signiﬁcant. The store parameter is negative from 7:00
to 11:00, while it is positive at 12:00, 15:00, 18:00, and 20:00. It means that while the parameteris negative during commuting and work hours, it is positive after work hours, indicating that theattractiveness of the facility changes with time. The park parameter is estimated to be signiﬁcantly
positive at 8, 16, 18, and 22 o’clock, indicating that parks improve travel comfort even duringcommuting hours. Therefore, it is essential to consider people’s behavioral choices that changewith time.

938.561205.75
1043.59
1331.941511.64
2013.911392.51
1568.692941.372941.92
1588.34
3825.68
1253.50
61.74
38.62

OD1
OD2
OD3
OD4
OD5
OD6
OD7
OD8
OD9
OD10
OD11
OD12
OD13
Initial KL
Final KL4669.465217.184736.165484.215383.706576.515641.30
5395.784201.684149.465461.615462.684977.52
62.11
39.02

-8.07
(-5.36)-0.61(-0.32)-0.02(-0.09)0.75(0.51)
0.36(1.09)-1.07(-1.03-7.26(-5.53)-1.02(-0.87)0.12(2.32)
0.39(0.72)
0.32(1.24)-1.23(-2.25)2488.102861.95
2537.363037.983642.564262.61
3121.913719.782965.362909.393752.225167.202903.40
59.35
38.47

9

8

7
estimation
value
Parameter(t-value)
Distance-7.27[θlength ](-5.81)Store-1.65[θstore ](-2.29)-0.03
Park(-0.29)[θpark ]Building
1.68[θbuilding ](1.69)
Station0.55[θstation ](3.08)-1.79
Road(-4.80)[θroad ]\time3427.723594.163481.723803.933614.234297.923790.804072.95
2792.79
2843.34
4114.315628.95
3136.36
64.34
39.50

-8.44
(-7.25)-1.56(-1.65)-0.03(-0.13)0.64(0.54)
0.46(2.49)-1.33(-2.65)10-7.69(-6.84)-1.35(-1.78)-0.01(-0.13)
1.21(1.23)
0.53(3.01)
0.53(-3.95)13-7.69(-6.84)-1.35(-1.78)-0.01(-0.13)
1.21(1.23)
0.53(3.01)-1.70(-3.95)14-8.21(-5.09)0.49(0.26)-0.31(-0.91)0.28(0.20)
0.13(0.40)-0.37(-0.36)15-7.79(-9.38)-1.56(-1.55)0.10(3.97)-0.03(-0.02)0.39(1.79)-1.08(-2.17)16-8.30(-5.52)1.50(1.86)
0.06(0.81)
0.08(0.09)-0.02(-0.12)0.11(0.22)

17-7.10(-5.63)
0.18
(-)
0.15
(6.32)
0.23(0.53)
0.12(1.13)
-0.51
(-)

18

3262.76 3613.72 2660.02 2660.02 7210.59 1601.23 3770.80 1073.25
3663.59 4184.69 2999.10 2999.09 7361.80 1705.23 4385.78 1274.89
3285.35 3535.54 2712.18 2712.18 7233.63 1673.55 3820.67 1501.34
3848.50 4478.87 3184.68 3184.68 7644.04 1776.02 4674.60 5625.85
3912.93 4857.01 4000.01 4000.01 7468.50 1611.32 4944.28 1262.44
5000.77 6122.09 4916.50 4916.50 8554.46 2033.07 5987.54 1589.62
4054.42 4790.06 3298.85 3298.85 8049.29 1752.36 5037.32 1379.25
3904.74 4840.61 4130.51 4130.51 7388.79 2251.69 4901.11 1274.53
2870.58 3786.83 3200.55 3200.55 6595.87 2658.06 4514.15 1168.57
2891.90 3801.75 3245.09 3245.09 6650.80 2711.00 4613.04 1168.26
3661.77 4908.80 4167.90 4167.90 7509.12 2299.69 4984.21 1280.00
6777.30 6707.38 6203.77 6203.77 7567.14 3679.07 6099.27 2820.88
3789.53 4801.31 3080.15 3080.15 6815.45 1387.24 4950.31 1246.91
67.00
63.94
64.62
64.63
63.77
65.69
61.82
62.25
38.50
39.11
40.46
40.46
42.31
40.52
38.35
41.77

-8.05
(-5.42)0.24(0.10)
0.05(0.34)
0.70(0.60)
0.22(0.49)-0.65(-0.43)12

TABLE 5 Estimation result at each hour

-9.15(-5.93)-0.07(-0.02)0.01(0.05)
0.68(0.62)
0.23(0.45)-0.60(-0.35)

112529.57
2904.222549.243033.303399.694380.103226.59
3901.232314.802328.823937.795483.872829.44
67.43
36.43

-9.56
(-6.35)-0.75(-0.81)-0.11(-0.43)0.70(0.63)
0.33(1.90)-0.99(-1.70)193543.863691.653601.843823.823681.844095.364088.763652.53
3337.143312.56
3721.396334.973759.31
62.50
37.84

-8.74
(-5.74)
1.64(1.27)
0.06(0.71)
0.38(0.30)-0.07(-0.28)0.44(0.57)

20-7.42(-5.77)-0.26(-0.55)0.13(4.83)
0.17(0.36)
0.39(2.32)-1.03(-2.11)

22-8.05(-5.19)0.29(0.14)-0.24(-0.72)0.14(0.10)
0.25(0.64)-0.72(-0.58)234815.87 1013.83 5232.84
5411.81 1107.167 5705.00
4819.12 1279.49 5252.22
5624.68 2695.44 5928.75
6513.47 1060.81 6005.75
7823.38 1367.89 7011.19
6116.29 1131.24 6193.52
6472.74 1095.94 5958.92
4819.28 1535.15 4961.88
4843.70 1755.09 5031.67
6524.03 1113.13 6009.58
6566.73 3240.76 6013.45
5354.96 1031.05 5198.24
61.89
63.46
65.63
38.96
40.25
42.86

-8.52
(-5.31)-0.09(-0.05)-0.27(-0.74)0.60(0.44)
0.20(0.55)-0.71(-0.63)21

Ogawa and Hato, and Ishii
16

Ogawa and Hato, and Ishii

17

SPATIAL SIMULATION USING SURROGATE MODELS
Using the estimated model, we attempt to plan the development of parks in the Toyosu area by considering the expected utility of tourists and congestion indices. We calculate the Pareto solution
set based on the park’s size and its location pattern, and conduct a Pareto frontier search. However, our RL model-based allocation calculation, which expresses the migration behavior, requires
recursive calculation using the Bellman equation, and the calculation cost is enormous. To solve
this problem, we construct a surrogate model using a feed-forward neural network to improve thespeed.
Design evaluation index
In this study, we consider the area of the park layer in the Toyosu area as a design variable. We
deﬁne the objective function of the design as the expected utility of a traveler per day as follows(24)
TS
1 23
(24)
∑ ∑ ∑ B(h, z j )S(t = 0; h, z j )Ball h=7z j ∈z t=1
Where B is the sum of OD matrixs for all periods, and B(h, z), is the OD matrix of OD pair z j
at time h, and S(t = 0; h, z j ) is the state value function of the start point of OD pair z j at timeh. We need to consider the above objective function maximization problem for all cells in thetarget regional network, but the computation of the optimal solution requires a large amount ofcomputational cost. Therefore, we search for a design variable that maximizes the index consideredby considering maintenance cost minimization as the objective function. The objective function of
the park improvement plan z( b, 1) is deﬁned as follows(25).

max .za,1 =min .zb,1 = ∑ (na − nmin
a )

(25)

a∈A

where na is the park area in the cell, and nmina is the current park area in the cell. We use the network
design problem as a multi-objective optimization of both equations(24 and 25) and try to improvethe computational efﬁciency in the next chapter.
Searching solution ﬂow
We show the searching solution ﬂow in Figure 8. First, we make the network in which the explanatory variables of the links are changed by random sampling and calculate the value of the objectivefunction. Next, we use the calculated objective function values and the corresponding network
data as the training data for the feed-forward neural network and train the model to speed up theallocation calculation to search for the Pareto solution. In this study, we obtain the number oftransitions between cells at each time and time step by performing allocation calculations based on
the parameters estimated for the test network. Then, we calculate the objective function and train
the model using a feed-forward neural network with the corresponding NW data as the train data.
(i)
The objective function (za,1 ) corresponding to each network G(i) is regarded as a Pareto solution
when the following conditions are satisﬁed(26).
(i)¬{∃G ∈ F, za,1 ≥ za,1 }where, F is the set of Pareto solutions.

(26)

Ogawa and Hato, and Ishii

18

FIGURE 8 Searching Pareto solution ﬂow using surrogate model

Training data and Model
This section creates a trainer that outputs the three objective function values described above basedon the allocation results for each period from 7:00 to 23:00 in the target site.
Training data
We used the estimated OD matrix and route choice models for each period to perform allocation
calculations on 6000 network data. The explanatory variables for each link are the estimated six,and the number of links is 272, so the input data is an array of 272 × 6. The output layer isrepresented by the expected utility of travelers per day z1 , the average trip length, and the Gini
coefﬁcient for population density considering time variation and the number of nodes in the output
layer.
NN ModelFirst, we show the structure of the NN used in this study in Figure 9. The input layer takes an arrayof 272 × 6 × 1 as input, and the output layer returns three outputs: expected utility, average trip
length, and Gini coefﬁcient for population density considering time variation. For NN learning, we
use the mean squared error (MSE) deﬁned between the prediction Xi and the correct answer X̂i as
the loss function. We set the learning rate to 0.001, the batch size to 20, and the number of epochs
to 20. The batch size is the number of data in one set when training data is divided, and the epochs
are the number of training sessions.

Ogawa and Hato, and Ishii

19

NW data

input: (272, 6, 1)

input: (272, 6, 256)

input: 208896

input: 128(272×6)

output: (272, 6, 256)

output: (272, 6, 128)

output: 128

output: 3

Dropout (0.5)

FlattenObjectivefunction(3)

Dropout (0.5)

FIGURE 9 Feed-forward neural network ﬂow

error/loss

ResultResult of speed up
Figure 10 shows the training process of the NN with the loss function value on the vertical axis andthe number of epochs on the horizontal axis. This ﬁgure shows the mean squared error and means
an absolute error in the validation check for each epoch and the mean absolute error for the training
data. It can be seen that the accuracy of the model improves with each epoch, and over-training
can be avoided.
We also compared the computation time of the objective function between the surrogate
model and the RL model on PC with Intel − Core(T M)i5 − 8250UCPU1.6GHz1.8GHz. The
computation time for one time period using the estimated OD matrix and spatial choice model for
100 network data was 440 seconds. On the other hand, the NN model required only 0.09 seconds,
which is 4888 times faster.

epoch

FIGURE 10 Learning process at each epoch of NN model

Result of planning park
Figure 11 shows the set of rejected NW and Pareto solutions. We can see the Pareto frontier, where
the expected utility increases with the increase in the park area. From the Pareto set in Figure 7.3,
we focus on Pareto solution A.

Ogawa and Hato, and Ishii

20

Figure 12 shows the relationship between park planning and OD points of Pareto solution
A. The Toyosu and Shin − Toyosu direction has the most extensive park area, but the OD matrix
from the T sukishima direction is signiﬁcant, so there is a need for park planning in the Harumi
area. It is thought that the development located close to the activity route of travelers is practical,
not only on the origin and destination points.

A
Pareto solution
Reject solution

FIGURE 11 Pareto frontier with minimum park area and maximum expected utility

FIGURE 12 Distribution of park planning area and OD matrix in Pareto solution A

Ogawa and Hato, and Ishii

21

CONCLUSION
We proposed a simultaneous estimation of OD matrix and route choice models, using link ﬂow orlink transition probabilities as observation and a surrogate modeling approach to speed up network
design.
In previous studies, estimation of the OD matrix and estimation of the parameters of thedemand model has been established as different methodologies. However, we have achieved integration of the estimation methods by making the OD matrix a latent variable.
Speciﬁcally, we deﬁned the simultaneous OD matrix and link ﬂow as a manifold and constructed a simultaneous estimation method by combining the em algorithm and NPL. This method
can guarantee the uniqueness of the solution from the viewpoint of information geometry. Furthermore, we compared the propose method with the conventional method and conﬁrmed the stabilityof the solution and the high accuracy of the estimation from numerical experiments.
And we attempted to solve the network design by constructing a surrogate model that replaces the high computational calculation cost. Using a surrogate model that learned input network pattern and output expected utility of the traveler are learned, we achieved a computational
speed 4888 times faster than the conventional framework. As a case study, we conducted a multiobjective optimization of the problem to minimize the area of park planning and maximize travelers’ expected utility. We conﬁrmed the Pareto frontier in the solution set with surrogate model.
In the future, it is necessary to perform the estimation for the whole period instead of anhourly estimation. It is a strong assumption that the values of the explanatory variables in the
network are independent and travel across time. And the fusion estimation method minimized thesum of the two KL divergences, but we cannot guarantee the uniqueness of the solution of thismethod. It is necessary to implement the estimation with guaranteed solution using real data fromdata handling. Due to the prolonged COVID-19, there has been a decrease in migration behavior
and OD matrix in the city, which has caused changes in the places where people are concentrated
and changes in behavior. It is necessary to clarify the extent of this decrease, and to develop urban
policy, urban design and evaluation for density control.
APPENDIX
Estimation method using link transition probability
We can use our simultaneous estimation method using link ﬂow and link transition probability.
Deﬁnition manifold when we observed link transition probability
First, we formulate the simultaneous distribution of OD trafﬁc and node transition probabilities
as a model manifold. When we transition to a node a ∈ A(k), the simultaneous distribution of its
transition probability and OD pair z j can be expressed as follows(27).
w j p(a|k, z j ; θ )p(k|z j ; θ )(27)p(k; θ )w j p(xl ; θ )=(28)∑ j ∑a w j p(xl ; θ )Next, we consider the observation manifold. We deal with the empirical distribution q̂(a|k)
of node transition probabilities as observations. In this case, the simultaneous distribution of OD
p(a, z j |k; xl) =Ogawa and Hato, and Ishii

22matrix and node transition probability can be expressed as follows(29).
q(a, z j |k) = q̂(a|k)q(z j |a, k)(29)
where, q(z j |a, k) is an arbitrary conditional distribution, deﬁned as the parameter corresponding toa, k, z j .
Estimation method using link transition probability
The KL divergence between the observed manifold and the model manifold is deﬁned as the objective function as follows(30).
q̂(a|k)q(z j |a, k)p(a, z j |k; ξ )k∈A a∈A z j ∈z
The parameter ξ = {θ , w} is estimated by KL divergence(30) minimization.

KL(q(z j |a, k)||ξ ) = ∑ ∑ ∑ q̂(a|k)q(z j |a, k) ln(30)θ̃ = min KL(q(z j |a, k)||ξ )

(31)

θ

w̃ = min KL(q(z j |a, k)||ξ )(32)w
Since this method uses transition probabilities as observations, it is impossible to convert OD trafﬁc
distribution to OD matrix. In this case, we need information about the ﬂow rate in the NW.
Next, we consider minimizing the KL divergence by the em algorithm.
e stepIn this step, we ﬁx Γ and consider the e-projection onto the observation manifold thatminimizes the KL divergence concerning q(z j |a, k). Since q(z j |a, k) is a probability distribution,
we can impose the constraint that ∑z j ∈z q(z j |a, k) = 1, and similar to the method when the transition
quantity is an observation, we can minimize the KL divergence q(z j |xl )(33).
q(z j |xl ) =p(z j |a, k; ξ )=p(z j |xl , ξ )p(xl , z j |ξ )=∑z j ∈z p(xl , z j |ξ )(33)

m step
Next, we ﬁx q(z j |a, k) and consider the m-projection onto the model manifold that minimizes the KL divergence with respect to ξ . This is equivalent to considering the following maximization problem(??).
L(ξ |q(z j |a, k)) = ∑ ∑ ∑ vl q(z j |a, k) ln p(a, z j |k; ξ )(34)k∈A a∈A z j ∈z

By using the proposed NPem algorithm, we can estimate the parameters that can guarantee theuniqueness of the solution.

Ogawa and Hato, and Ishii

23

Fusion estimation using link ﬂow and transition probability
In this study, we use the described PP data and mobile spatial statistics data. From the PP data,
we can obtain the transition probability distribution between cells. In addition, we can obtainthe number of cell states in 500m mesh units from the mobile spatial statistics data. Thus, we
can obtain the number of transitions between cells in a 500m mesh NW. In this study, the multiscale network handle meshes of 100m and 50m, so we observe the cell transition probabilities
as observations for scales of 100m or less. Therefore, we will explain the estimation of fusion
estimation method that uses the transition amount and transition probability as observations.
We deﬁne the objective function as the sum of the KL divergence of the link ﬂow as anobservation and the link transition probability as an observation(35).
min KL(q(z|x)||ξ ) + KL(q(z|a, k)||ξ )(35)

θ ,q(z|a,k)This can be interpreted as a problem of deﬁning a model manifold with the same coordinate system and an observation manifold with the same coordinate system, and minimizing the degree ofseparation between each manifold, from q(z j Γa, k) = q(z j |xl ).
In the following, each step of the em algorithm for minimizing the KL divergence is explained.
e stepFix ξ and minimize the KL divergence with respect to q(z j Γa, k). We can use the Lagrangian undecided multiplier k,a and express it as follows(36).
L(q(z j |a, k), λ |ξ ) = ∑ ∑ ∑ q̂(a|k)q(z j |a, k) ln
k∈A a∈A z j ∈zq̂(a|k)q(z j |a, k)p(a, z j |k; ξ )q̂(xl )q(z j |xl )− ∑ λk,a ( ∑ q(z j |a, k) − 1)p(xl , z j |ξ )z j ∈zxl ∈E z j ∈zk,a∈A
From KKT conditions and q(z j |k, a) = q(z j |xl ),+ ∑ ∑ q̂(xl )q(z j |xl ) ln

(36)q(z j |a, k)∂ L(q(z|x), λ |ξ )= {q̂(a|k) + q̂(xl )} ln+ ηk,a(37)∂ q(z j |xl )p(z j |a, k; ξ )where ηk,a is a constant. From the above, when ξ is ﬁxed, the q(z j Γa, k) that minimizes the KL
divergence is the following equation.
m step
Next, we ﬁx q(z j Γa, k) and consider the m-projection onto the model manifold that minimizes the KL divergence with respect to ξ . It is equivalent to considering the following maximization problem.
L(ξ ||q(z j |a, k)) = ∑∑ q̂(a|k)q(z j |xl ) ln p(a, z j |k; ξ ) + q̂(xl )q(z j |xl ) ln p(xl , z j |ξ )(38)k,a∈A z j ∈z

Note that this maximization problem is not a single KL divergence, but a projection problem dealing with the sum of two KL divergences, so it is difﬁcult to discuss ﬂatness.

Ogawa and Hato, and Ishii

24

REFERENCES
1.
Sakai, H., M. Shimizu, T. Yoshimura, and E. Hato, Psychological Reactance to Mobility
Restrictions Due to the COVID-19 Pandemic: A Japanese Population Study. Frontiers in
Psychology, 2021.
2.
Hato, E., Development of behavioral context addressable loggers in the shell for travelactivity analysis. Transportation Research Part C: Emerging Technologies, Vol. 18, 2010,
pp. 55–67.
3.
Oyama, Y. and E. Hato, Link-based measurement model to estimate route choice parameters in urban pedestrian networks. Transportation research part C: emerging technologies,
Vol. 93, 2018, p. 62–78.
4.
Liu and Fricker, Estimation of a trip table and the θ parameter in a stochastic network.
Transportation Research Part A: Policy and Practice, Vol. 30, 1996, p. 287–305.
5.
Dial and Robert, A probabilistic multipath trafﬁc assignment model which obviates pathenumeration. Transportation research, Vol. 5, No. 2, 1971, p. 83–111.
6.
Yang, Meng, and M. GH, Simultaneous estimation of the origin-destination matrices andtravel-cost coefﬁcient for congested networks in a stochastic user equilibrium. Transportation science, Vol. 35, 2001, p. 107–123.
7.
Balakrishna, Ben-Akiva, and Koutsopoulos, Ofﬂine calibration of dynamic trafﬁc assignment: simultaneous demand-and-supply estimation. Transportation Research Record, Vol.
2003, 2007, p. 50–58.
8.
Wang, Ma, Yong, K. Gong, K. C. Henricakson, M. Xu, and Y. Wang, Ofﬂine calibrationof dynamic trafﬁc assignment: simultaneous demand-and-supply estimation. PloS one,Vol. 11, 2016.
9.
Lo and Chan, Simultaneous estimation of an origin–destination matrix and link choiceproportions using trafﬁc counts. Transportation Research Part A: Policy and Practice,Vol. 37, 2003, p. 771–788.
10.
Oyama, Y. and E. Hato, Prism-based path set restriction for solving markovian trafﬁc
assignment problem. Transportation Research Part B: Methodological, Vol. 122, 2019, p.
528–546.
11.
Fosgerau, Mogens, Frejinger, Emma, Karlstrom, and Anders, A link based network route
choice model with unrestricted choice set. Transportation Research Part B: Methodological, Vol. 56, 2013, pp. 70–80.
12.
Bhat, C. R., An endogenous segmentation mode choice model with an application to intercity travel. Transportation science, Vol. 31, 1997, p. 34–48.
13.
Vij, A., A. Carrel, and J. L. Walker, Incorporating the inﬂuence of latent modal pref- erences on travel mode choice behavior. Transportation Research Part A: Policy and Practice, Vol. 54, 2013, p. 164–178.
14.
Akaho, Information Geometry in Machine Learning. Journal of the Society of Instrument
and Control Engineers, Vol. 44, No. 5, 2005, pp. 299–306.
15.
Csiszaér and G. Tusnaédy, Information geometry and alternating minimization problems.
Statistics Decision, Supplement Issue, , No. 1, 1984, pp. 205–237.
16.
Amari, Information geometry of the EM and em algorithms for neural networks. Neuralnetworks, Vol. 8, No. 9, 1995, pp. 1379–1408.
17.
Fosgerau, Mogens, Melo, Emerson, Palma, A. de, Shum, and Matthew, Discrete choice
and rational inattention: A general equivalence result. Available at SSRN 2889048, 2017.

Ogawa and Hato, and Ishii
18.

19.

20.

21.
22.

23.

24.
25.

25

Thelin, J. Salmon, S. Gorrell, S. Bunnell, G. Bird, C. Ruoti, E. Selin, and J. Calogero, singsurrogate models to predict nodal results for fatigue risk analysis. International Journal of
Fatigue, Vol. 146, 2021.
Mai, J. Kang, and J. Lee, A machine learning-based surrogate model for optimization oftruss structures with geometrically nonlinear behavior. Finite Elements in Analysis andDesign, Vol. 196, 2021.
Chang, K.-H. and C.-Y. Cheng, Learning to Simulate and Design for Structural Engineering. Proceedings of the 37 the International Conference on Machine Learning,PMLR, Vol.
119, 2020.
Hao, V. G. . J.-K., Transit network design and scheduling: A global review. Transportation
Research Part A: Policy and Practice, Vol. 42, 2008, p. 1251–1273.
Farahani, R. Z., E. Miandoabchi, W.Y.Szeto, and H. Rashidid, A review of urban transportation network design problems. European Journal of Operational Research, Vol. 229,
2013, p. 281–302.
Oyama, Y. and E. Hato, Pedestrian Activity-based Network Design based on MultiObjective Programming. Journal of City Planning of institute Japan, Vol. 52, No. 3, 2017,
pp. 810–817.
Vliet, Selected node-pair analysis in dial s assignment algorithm. Transportation ResearchPart B: Methodological, Vol. 15, 1981, p. 65–68.
Aguirregabiria, V. and P. Mira, Swapping the Nested Fixed Point Algorithm: A Classof Estimators for Discrete Markov Decision Models. The Econometric Society, Vol. 70,
No. 4, 2002, pp. 1519–1543.