1 Introduction

Refugees have historically always been part of human migration, seeking locations of greater safety and security for themselves and their families. Now, with the increasing impacts of climate change, and associated increases in the severity, as well as the frequency, of many disasters, including hurricanes and floods, coupled with numerous conflicts, strife, and violence in many parts of the globe, the number of refugees and asylum seekers is escalating. According to the United Nations [52], the number of international migrants was approximately 258 million (3.4% of the global population) in the year 2017, with the total number of international migrants increasing by almost 50% since 2000 (United Nations [52]). Also, since the new millennium, the number of refugees and asylum seekers has increased from 16 to 26 million, constituting approximately 10% of the international migrants. And, according to the United Nations Refugee Agency [55], the global population of forcibly displaced individuals increased by 2.3 million people in 2018. Almost 70.8 million individuals had to move due to persecution, conflict, violence, or human rights violations. Now with the COVID-19 pandemic, declared by the World Health Organization on March 11, 2020, Dolmans et al. [12] report that COVID-19 is likely to exacerbate what is already a humanitarian emergency in terms of a global refugee crisis. The authors argue that refugees may encounter increased difficulty in seeking asylum due to measures imposed by governments in response to the pandemic.

Jackson [16] reports that the legal definition of a refugee was provided in 1951 by the United Nations Convention Relating to the Status of Refugees. The Convention defined a refugee as an individual living outside his or her country of nationality, who is unable or unwilling to return because of a well-substantiated fear of persecution due to race, religion, nationality, membership in a political social group, etc. Here we also consider humans adversely affected by climate change, as refugees, and note that, as emphasized by Hebert et al. [13], among the most studied causes of human migration are climate issues and conflicts, as well as economic reasons. In this paper, since the focus is on refugees, we do not emphasize the latter, but concentrate on the two former reasons for migration.

The news media in recent years has been filled with reports of refugees traveling from points of origin towards destinations that they ultimately hope to reach, sometimes using challenging and even dangerous modes of transportation. For example, as emphasized in Nagurney and Daniele [31], the United Nations Refugee Agency [54] reported a maritime refugee crisis with 137,000 crossing the Mediterranean Sea to Europe, in the first half of 2015, using very risky transport modes, and with many others unsuccessfully attempting such a crossing. 800 lost their lives in the largest refugee shipwreck on record that spring.

The United Nations International Organization for Migration notes an article by Conant [5] in the National Geographic highlighting, through vivid graphical maps, several of the major mass migration routes that were then underway. The article emphasizes that those fleeing the Middle East and Africa to Europe are not the only ones traveling along evolving and dangerous migration routes. Many other countries are also being characterized by an exodus of their citizenry and/or are struggling with becoming transit points, smugglers’ routes, or desired destination nodes for migrants. The routes, which may be secretive, are cutting paths through Central America and Mexico; the Horn of Africa (with almost one million Somali refugees), and other countries including: Bangladesh, Myanmar, and Malaysia, along with the headline-grabbing migrations through East Africa and the Mediterranean Sea.

Governments of many nations on multiple continents are faced with the challenges associated with the migration of refugees. Various strategies are being utilized to regulate the flow to various destinations. Helbling and Leblang [14] eruditely emphasize that migration regulations are central to a nation’s sovereignty since they directly impact the makeup of the community in a country. The United Nations [51], in turn, reports that migration policies in both origin and destination countries play a crucial role in determining the migratory flows. We can expect refugee migratory flows to continue to increase, posing a critical need for the provision of rigorous tools for policy-makers and decision-makers for the quantification of refugee migratory flows and the impacts of various regulations. That is one of the goals of this paper.

Indeed, the development of computable mathematical models that can assist in the prediction of refugee flows, along with the assessment of various regulations, are sorely needed. Interestingly, although there has been a growing literature on network models of human migration, certain features, alluded to above, have been missing. We first recall the relevant literature in Sect. 2, followed by a delineation of the specific, concrete contributions in this paper.

2 Literature review of human migration network models and our contributions

Multiple disciplinary approaches, as well as perspectives, have been applied to the study of human migration, including sociological, political, economic, as well as operations research based, with associated methodologies. Please refer to Rahmati and Tularam [45] for a review of a spectrum of frameworks for the analysis of human migration models. Recently, there has been a renewed interest in the development of network models for human migration, with an example being the work of Nagurney and Daniele [31], which is undoubtedly the first to include regulations in the modeling of international human migration networks. That work builds on the earlier contributions of Nagurney [24], who developed a multiclass migration network model but without migration costs, and that of Nagurney [25] and Nagurney, Pan, and Zhao [35, 36], who introduced flow-dependent migration costs, with the latter work capturing class transformations.

Pan and Nagurney [41] proposed a multi-stage, single class, Markov chain model in order to model chain migration, which was not a feature in the above-noted research, and provided theoretical results connecting a sequence of variational inequalities and a non-homogeneous Markov chain. Later, Pan and Nagurney [42] applied the theory of evolution variational inequalities to model the dynamic adjustment processes associated with human migration. For some background on finite-dimensional variational inequalities, we refer the interested reader to the book by Nagurney [26], and for background on evolution variational inequalities, see the book by Daniele [11]. In the network models of human migration, migrants reflect their preferences for destination versus origin nodes through utility functions, which are distinct for each class of migrant, and the utility functions depend on the populations at the different locations. Other extensions have included the inclusion of uncertainty in the utility and migration cost functions, as well as the initial populations, as in the work of Causa, Jadamba, and Raciti [4]. Kalashnikov et al. [18] proposed a human migration model with a conjectural variations equilibrium (CVE), which extended the framework of Isac, Bulavsky, and Kalashnikov [15] in that the conjectural variations coefficients are continuously differentiable functions of the total population at the destination and of the group’s fraction in it. Cappello and Daniele [3], subsequently, constructed a Nash equilibrium model of human migration with features of conjectural variations. The utility and migration cost functions for a given class depend on the populations and flows, respectively, of that class. They also provided solutions numerical examples inspired by the refugees flows from Africa through the Mediterranean Sea to Italy in the year 2018.

Nevertheless, none of the earlier models noted above include regulations, except for that of Nagurney and Daniele [31]. Furthermore, all the above models assumed that a path from an origin node to a destination node consists of a single link, whereas, in reality (cf. Conant [5]), migrants, in the form of refugees, in particular, may actually be selecting paths, consisting of multiple links. Moreover, different pathways even within an intermediate country may be used to their ultimate destination. Some locations, in turn, may act as meeting points for refugees, signifying that there are, typically, different paths/routes that a refugee may take from an origin node to a destination node. Interestingly, although having multiple path options from origins to destinations has been typical in transportation network models (cf. Beckmann, McGuire, and Winsten [2], Dafermos [7, 8], Sheffi [47], Zhang and Nagurney [56], among many others), as well as in spatial price equilibrium models, in which the flows are commodities (see Dafermos and Nagurney [9], Nagurney and Aronson [29]) and in supply chain models (cf. Nagurney [27, 28]), Nagurney, Yu, and Besik [37], Nagurney et al. [38], Toyasaki, Daniele, and Wakolbinger [50], and the references therein), including multiple routes has not been done in the context of human migration network models. We fill this research gap in this paper. By allowing for this important feature we then capture additional possible regulations via novel constraints.

We also establish that the model without regulations is isomorphic to the well-known traffic network equilibrium (TNE) problem with fixed demands (cf. Beckmann, McGuire, and Winsten [2], Dafermos [7], Dafermos and Sparrow [10], Sheffi [47], Nagurney [26], Patriksson [43]) through a supernetwork (see Nagurney and Dong [32]) equivalence. This identification, which also advances the theory of network equilibria as a unifying mechanism for many problems, enables the transferral of numerous algorithms that have been developed for TNE to general migration network problems. We note that multitiered financial equilibrium problems had been transformed into fixed demand traffic network equilibrium problems (Liu and Nagurney [21]); supply chain network equilibrium problems have been as well (Nagurney [27]), as have electric power generation and distribution network problems (Nagurney et al. [34]). However, in this paper, the first such transformation is done for a migration equilibrium problem. Moreover, the refugee migration network problem is defined over a general network topology. The refugee migration model with regulations can, hence, be transformed into a TNE with side constraints (see [17]; Larsson and Patriksson [20]).

The specific contributions to the literature of this paper are:

  1. 1.

    An international human migration network model is constructed, which allows for route choices by the migrants, which are refugees.

  2. 2.

    The routes consist of one or more links, with cost functions that capture congestion, a factor that has been seen in practice.

  3. 3.

    The model is then extended to include regulations that can be imposed by distinct multiple countries. In previous work (cf. Nagurney and Daniele [31]), it was assumed that a single country imposes the regulations on migrants.

  4. 4.

    A supernetwork transformation into a traffic network equilibrium problem with fixed demands is constructed. This identification enables the transfer of algorithmic schemes for the TNE problem, which has had a long history, to the novel application domain of refugee/migration networks.

  5. 5.

    Theoretical results are presented, along with an algorithm, with nice features for computations.

  6. 6.

    Numerical examples illustrate the modeling and algorithmic framework and enable insights for policy-makers and decision-makers.

The paper is organized as follows. In Sect. 3, the multiclass, multipath refugee models without and with regulations are developed. In Sect. 4, the supernetwork transformation of the problem into a traffic network equilibrium problem is presented. In Sect. 5, we propose an algorithm, which is then applied to a series of numerical examples. In Sect. 6, we summarize our results, present our conclusions, and provide suggestions for future research.

3 The multiclass, multipath refugee migration models

In this section, we construct refugee migration models consisting of multiple classes of refugees and with the possibility of multiple routes/paths between an origin node and a destination node associated with refugee flows. We first present the model without regulations and then the one with. The regulations can be imposed by multiple governments.

The refugee migration network consists of a graph \(G=[N,L]\), where N denotes the set of nodes and L the set of directed links representing the direction of flow. Origin nodes correspond to initial locations of the refugees and there could be multiple such nodes in a country. Similarly, destination nodes correspond to destinations in countries that the refugees may be interested in traveling to. We also consider the situation that the refugees may elect to stay in their own country for various reasons, including that the cost of migration may exceed the net utility gain. A typical origin node is denoted by i and a typical destination node by j, with the origin/destination (O/D) pair of nodes associated with origin i and destination j denoted by \(w_{ij}\). The set of routes/paths joining O/D pair \(w_{ij}\) is denoted by \(P_{w_{ij}}\), with a typical route referred to by r. The set of all routes in the refugee network is denoted by P where \(P=\cup _{i,j=1}^n P_{w_{ij}}\). There are J classes of refugees, with a typical class denoted by k. According to the United Nations High Commissioner for Refugees [53], a refugee has a justifiable fear of persecution because of race, religion, nationality, political opinion, or belonging to a specific social group. Depending on the specific scenario and refugee application, one can define a class accordingly and even consider other identifiable attributes such as the social-economic level, professional status, etc. Here we consider a class as being distinct in terms of associated utility at origin and destination nodes as well as cost associated with migration. In creating the model, for the sake of flexibility, we do not restrict ourselves to specific attributes. We see that an O/D pair of nodes corresponds to nodes in distinct countries, except in the case of a pair of nodes in the same country corresponding to the same location and with a route consisting of a single link joining such a pair of nodes, signifying the decision to stay (and not migrate). We assume, for simplicity, that there are n origin nodes and n possible destination nodes. In our models we assume that there is no repeat, that is, chain, migration. Also, all vectors are column vectors.

Table 1 Common notation for the refugee migration models

The common notation for the models is given in Table 1. Observe that, according to Table 1, the utility functions of each refugee class can, in general, depend upon the population distribution of all the refugees in all the locations; moreover, the cost on a route/path can, in general, depend upon the vector of all the refugee flows. Of course, the specific structure of these functions would be adapted to the migration scenario(s) under study.

We first present the conservation of flow equations. Since the route flows must be nonnegative, we have that

$$\begin{aligned} x_r^k\ge 0, \quad \forall r\in P, \forall k. \end{aligned}$$
(1)

Furthermore, the refugee flows out of an origin node i must satisfy:

$$\begin{aligned} {\bar{p}}_i^k=\sum _{j=1}^n\sum _{r\in P_{w_{ij}}} x_r^k,\quad \forall i,\forall k. \end{aligned}$$
(2)

Note that Eq. (2) states that for each class of refugee at each origin node, the sum of the flows out by the class must be equal to the total initial population of that class at the origin node. Remember that refugees can also elect to stay and that would correspond to O/D pairs \(w_{ii}\); \(i=1,\ldots ,n\), with the corresponding flows on the paths consisting of single links associated with such O/D pairs.

Furthermore, the volume of population of each class k at each destination node j, after migration takes place, must satisfy the following equation:

$$\begin{aligned} p_j^k=\sum _{i=1}^n\sum _{r\in P_{w_{ij}}}x_r^k,\quad \forall j,\forall k. \end{aligned}$$
(3)

Moreover, we have the following conservation of flow equations which relate the link flows to the route flows:

$$\begin{aligned} f_a^k=\sum _{r\in P} x_r^k\delta _{ar}, \quad \forall a,\forall k, \end{aligned}$$
(4)

where \(\delta _{ar}=1\), if link a is contained in route r and 0, otherwise. According to (4), the flow on a link of a class is equal to the sum of flows of that class on routes that use the link.

In view of the conservation of flow equations (4), we may define link cost functions in route/path flows, such that \({\hat{c}}_a^k={\hat{c}}^k_a(x)\equiv c_a^k(f)\), for all links a and for all classes of refugees k.

The cost on a route r is equal to the sum of costs on the links that make up the route, that is,

$$\begin{aligned} C_r^k(x)=\sum _{a\in L} {\hat{c}}_a^k(x)\delta _{ar},\quad \forall k,\forall r. \end{aligned}$$
(5)

We define the feasible set \(K^1\equiv \{(p,x)\vert x\ge 0,\,\text{ and }\, (2)\,\text{ and }\, (3)\,\text{ hold }\}\).

A sample refugee network is depicted in Fig. 1. Note that the direct link from an origin node to a destination node with the same numeric label will have, once the model is solved, the volume of the population that remains at that location (and does not migrate). The costs on such links are set equal to zero, since those links are associated with staying at the origin location(s).

Fig. 1
figure 1

Sample refugee network topology

3.1 Equilibrium conditions for the multiclass, multipath refugee migration model without regulations

We first present the equilibrium conditions for the multiclass, multipath refugee migration model without regulations, accompanied by variational inequality formulations. We then introduce its counterpart with regulations. The equilibrium conditions are generalizations of classical migration equilibrium conditions of Nagurney, Pan, and Zhao [35] to include multiple possible routes.

Definition 1

(Multiclass, Multipath Refugee Migration Equilibrium without Regulations)

A vector of populations and refugee migration flows \((p^*,x^*)\in K^1\) is in equilibrium if it satisfies the following conditions: For each class k; \(k=1, \ldots ,J\), and each pair of origin/destination nodes ij; \(i,j=1, \ldots ,n\), and all routes \(r\in P_{w_{ij}}\) we have that

$$\begin{aligned} u_i^k(p^*)+C_r^k(x^*) {\left\{ \begin{array}{ll} =u_j^k(p^*)-\lambda _i^{k*}, &{} \text{ if }\, x_{r}^{k*}>0,\\ \ge u_j^k(p^*)-\lambda _i^{k*}, &{} \text{ if } \, x_{r}^{k*}=0,\end{array}\right. } \end{aligned}$$
(6)

and

$$\begin{aligned} \lambda _i^{k*}{\left\{ \begin{array}{ll} \ge 0, &{}\text{ if } \, \sum _{j=1}^n\sum _{r\in P_{w_{ij}}, j\ne i}x_r^{k*}=\bar{p}_i^k,\\ = 0, &{} \text{ if } \, \sum _{j=1}^n\sum _{r\in P_{w_{ij}}, j\ne i}x_r^{k*} <{\bar{p}}_i^k. \end{array}\right. } \end{aligned}$$
(7)

We assume, as noted earlier, that \(C_r^k=0\), for all r and k, if r consists of a route with an origin node and a destination node in the same country/location, since those who decide to stay in their origin countries and not migrate encumber a migration (travel) cost of zero.

Hence, according to (6), refugees will select routes from their origin nodes to their destinations such that the net gain in utilities for each class is maximal. Also, as emphasized in Nagurney and Daniele [31], but for international human migration network models in the case of a single path/link between each pair of nodes, \(\lambda _i^{k*}\) is necessary since the level of the population \(\bar{p}_i^k\) may not be sufficiently large so that the gain in utility \(u_j^k-u_i^k\) is exactly equal to the migration cost. Nevertheless, the utility gain minus the migration cost will be maximal and nonnegative, and the net gain will be equalized for all countries and classes which have a positive flow out of a node in a country (see also Nagurney, Pan, and Zhao [35]). For equilibrium conditions for spatial price equilibrium problems and transportation network equilibrium problems, see Beckmann, McGuire, and Winsten [2], Dafermos and Sparrow [10], Daniele [11], Nagurney [26], Nagurney, Besik, and Dong [30], Nagurney, Li, and Nagurney [33], Samuelson [46], Takayama and Judge [49], Zhang and Nagurney [56], and the references therein.

The variational inequality formulation of the multiclass, multipath refugee migration equilibrium conditions according to Definition 1 is given in Theorem 1 below.

Theorem 1

(Variational Inequality Formulation of the Refugee Migration Model without Regulations in Path Flows)

A population and refugee flow pattern \((p^*,x^*)\in K^1\) is a refugee migration equilibrium without regulations according to Definition 1, if and only if it satisfies the variational inequality problem in path flows

$$\begin{aligned} -\langle u(p^*), p-p^*\rangle + \langle C(x^*), x-x^*\rangle \ge 0, \quad \forall (p,x) \in K^1, \end{aligned}$$
(8)

where \(\langle \cdot ,\cdot \rangle \) denotes the inner product in the appropriately dimensioned Euclidean space.

Proof

See the Appendix.

Existence of at least one solution to variational inequality (8) follows from the standard theory of variational inequalities (see Kinderlehrer and Stampacchia [19] Theorem 3.1) under the assumption of continuity of the utility functions u and the migration cost functions c, since the feasible convex set \(K^1\) is compact.

Alternative variational inequality formulations can induce distinct algorithmic schemes. Hence, for completeness, we now provide a link flow variational inequality formulation equivalent to the path flow one in (8). We first define the feasible set \(K^2\equiv \{(p,f) \vert \exists x \,\text{ such } \text{ that }\, (1){-}(4)\, \text{ hold }\}\). The proof is immediate by making use of the conservation of flow equations (4), which relate the link flows to the path flows. \(\square \)

Corollary 1

(Variational Inequality Formulation of the Refugee Migration Model without Regulations in Link Flows)

A population and refugee link flow pattern \((p^*,f^*)\in K^2\) is a refugee migration equilibrium without regulations according to Definition 1, if and only if it satisfies the variational inequality problem in link flows

$$\begin{aligned} -\langle u(p^*), p-p^*\rangle + \langle c(f^*), f-f^*\rangle \ge 0, \quad \forall (p,f) \in K^2. \end{aligned}$$
(9)

We note that although alternative variational inequality formulations of both traffic network equilibrium problems as well as spatial price equilibrium problems have been constructed (cf. Nagurney [26] and the references therein), the above are the first such constructs for problems of human migration with multiple possible paths between O/D pairs.

3.2 Variational inequality formulation of the refugee migration model with regulations

We now construct general constraints of relevance to regulation of refugee migratory flows. Unlike in the model of Nagurney and Daniele [31], where it was assumed that only a single country was imposing regulations, here we allow for multiple countries to impose regulations. In addition, the set of constraints is of relevance to refugees in transit through a country that is neither an origin node for them nor a destination node.

We denote a specific country by h, where \(h=1,\ldots ,H\). Also, we define the set \(O^h\) consisting of origin nodes of refugees from countries/locations that the country h imposes a regulation on, and we let \(D^h\) denote the set of destination nodes, which lie in country h. Also, \(C^h\) denotes the set of refugee classes that country h imposes the regulations on and \(U^h\) is the nonnegative upper bound imposed by country h on refugee migratory flows.

The constraints can then be stated as follows:

$$\begin{aligned} \sum _{i\in O^h} \sum _{j\in D^h}\sum _{k\in C^h}\sum _{r\in P_{w_{ij}}}x_r^k\le U^h, \quad h=1,\ldots ,H. \end{aligned}$$
(10)

Observe that the set of constraints (10) is sufficiently general to capture specific, distinct migration regulations in practice. such as: an upper bound (which may be zero) of all classes from a certain country or countries; or an upper bound on a single class or several classes from a specific country or countries.

Furthermore—and this attests to the flexibility and generality of our multiclass, multipath framework—in the case that there are restrictions on the migratory flows on transit links, one can identify all such paths that utilize such links and restrict the total flow on them accordingly.

For the refugee migration model with regulations, the equilibrium conditions (6) and (7) are still relevant but with a new feasible set \(K^3\) defined as below to include the constraints (10):

$$\begin{aligned} K^3\equiv K^1\cap \{x\vert (10) \,\text{ is } \text{ satisfied }\}. \end{aligned}$$
(11)

When the feasible set \(K^3\) is nonempty then it follows immediately (cf. Nagurney [26]) that the variational inequality for the refugee migration model with regulations in path flows is as given in Theorem 2 below. We remark that, since we do not expect regulations to be imposed on refugees in the origin node country by the origin node country, refugees can remain in the country even if all other pathways are unavailable to them and, hence, \(K^3\) is nonempty.

Theorem 2

(Variational Inequality Formulation of the Refugee Migration Model with Regulations in Path Flows) A population and refugee migration flow pattern \((p^*,x^*)\in K^3\) is a refugee migration equilibrium with regulations in path flows, if and only if it satisfies the variational inequality problem

$$\begin{aligned} -\langle u(p^*), p-p^*\rangle + \langle C(x^*), x-x^*\rangle \ge 0, \quad \forall (p,x) \in K^3. \end{aligned}$$
(12)

Note that variational inequality (12) differs from variational inequality (8) in that the feasible set \(K^3\) is distinct from \(K^1\).

One can also construct the analogue to variational inequality (9), which is in link flows, but with regulations, by including the regulatory constraints into the expanded feasible set.

We now provide a variational inequality equivalent to the one in (12) but in path flows only (and not in the population distribution and path flows). We define a new feasible set \(K^4\equiv \{x\vert x\in R_+^{Jn_P}, \text{ and }\, (2),\,\text{ and }\, (1)\, \text{ hold }\}\). Because of the form of the utility functions as defined in Table 1, and using (3), we can define utility functions \({\hat{u}}_i^k(x)\), \(\forall i\), \(\forall k\), such that

$$\begin{aligned} {\hat{u}}_i^k(x)\equiv u_i^k(p). \end{aligned}$$
(13)

Note that variational inequality (12) can now be expressed solely in path flows. Indeed, the first term in (12) is:

$$\begin{aligned}&-\sum _{k}\sum _{i}u_i^k(p^*)\times (p_i^k-p_i^{k*}) \nonumber \\&\quad =-\sum _{k}\sum _{i}{\hat{u}}_i^k(x^*)\times \left( \sum _l\sum _{r\in P_{w_{li}}}x_r^k-\sum _l\sum _{r\in P_{w_{li}}}x_r^{k*}\right) . \end{aligned}$$
(14)

Hence, it follows that variational inequality (12) is equivalent to the variational inequality: determine \(x^*\in K^4\) such that

$$\begin{aligned} \sum _k\sum _{i}\sum _{j}\sum _{r\in P_{w_{ij}}}(-{\hat{u}}_j^k(x^*)+ C_{r}^k(x^*))\times (x_{r}^k-x_{r}^{k*})\ge 0,\quad \forall x\in K^4. \end{aligned}$$
(15)

We now put variational inequality (15) into standard form (cf. Nagurney [26]): determine \(X^*\in {{{\mathcal {K}}}}\) such that

$$\begin{aligned} \langle F(X^*), X-X^*\rangle \ge 0,\quad \forall X\in {{{\mathcal {K}}}}, \end{aligned}$$
(16)

where F is a given continuous function from \({{{\mathcal {K}}}}\) to \(R^N\), \({{{\mathcal {K}}}}\) is a given closed convex set.

Indeed, we can set \({{{\mathcal {K}}}}\equiv K^4\), \(X\equiv x\), and \(N=n_P\). Also, we can define the vector F as consisting of the elements: \(-{\hat{u}}_j^k(x)+C_r^k(x)\), for all origin nodes i, for all destination nodes j, for all routes \(r\in P_{w_{ij}}\), and for all classes k, and the conclusion follows.

3.3 Illustrative examples

We now proceed to illustrate the above results in several simple examples. The refugee network consists of the topology depicted in Fig. 2.

There are two origin country nodes and the same two country destination nodes. According to the network in Fig. 1, refugees residing in country 1 are not interested in migrating to country 2. On the other hand, refugees residing in country 2 are interested in the possibility of migrating to country 1. There are two available paths joining country 2 with country 1. The links are labeled in Fig. 1.

Fig. 2
figure 2

Network topology for illustrative refugee migration examples

The routes are comprised of links and are enumerated as follows:

$$\begin{aligned} r_1=(1), \quad r_2=(2),\quad r_3=(3,4), \quad r_4=(5,6,7). \end{aligned}$$

We consider a single class of migrant. Consequently, we suppress the superscript 1 in the notation. The data are as follows: \({\bar{p}}_1=100\) and \({\bar{p}}_2=200\) with the utility functions: \(u_1(p)=-p_1+1000\) and \(u_2(p)=-p_2+500\). Note that, according to these utility functions, country 1 is more attractive to the refugees than country 2, under no population of refugees, because of the term of 1000 in the utility function associated with country 1 as opposed to the fixed term of 500 in the utility function associated with country 2.

The link migration cost functions are:

$$\begin{aligned} c_1= & {} c_2=0, \\ c_3(f)= & {} f_3+200, \quad c_4(f)=f_4+100, \\ c_5(f)= & {} f_5+30,\quad c_6(f)=.5f_6+40,\quad c_7(f)=f_7+30. \end{aligned}$$

We see that route \(r_2\) consists of two legs of a journey, whereas route \(r_3\) is comprised on three legs. However, the fixed components of the migration cost on the legs of the former are much higher than those on the latter.

We first consider the case without regulations. It is easy to compute the equilibrium solution, using simple algebra. Indeed, we find that:

$$\begin{aligned} x_{r_1}^*=100,\quad x_{r_2}^*=75, \quad x_{r_3}^*=25, \quad x_{r_4}^*=100; \end{aligned}$$

hence,

$$\begin{aligned} p_1^*=225,\quad p_2^*=75, \end{aligned}$$

with associated utilities being:

$$\begin{aligned} {\hat{u}}_1(x^*)=u_1(p^*)=775,\quad {\hat{u}}_2(x^*)=u_2(p^*)=425, \end{aligned}$$

and the incurred migration costs on the routes at equilibrium: \(C_{r_1}=0\), \(C_{r_2}=0\), \(C_{r_3}=C_{r_4}=350\). Moreover, \(\lambda _1^*=\lambda _2^*=0\).

It is clear that the equilibrium conditions (6) and (7) hold.

Observe that, initially, before the refugee migration takes place, \(u({\bar{p}}_1)=900\) and \(u_2({\bar{p}}_2)=300\). Once equilibrium is achieved, those who have migrated from country 2 to country 1 more than double their utility, whereas those who remain in country 2 experience a gain in utility of over 33%. Those in country 1, because of the increase in the number of refugees and that they are not migrating, suffer a reduction in utility of approximately 14%.

We now suppose that a regulation is imposed on destination node 1 by country 1 of the following form:

$$\begin{aligned} x_{r_3}+x_{r_4}\le U^1=25. \end{aligned}$$

The new equilibrium solution is:

$$\begin{aligned} x_{r_1}^*=100,\quad x_{r_2}^*=175, \quad x_{r_3}^*=0, \quad x_{r_4}^*=25; \end{aligned}$$

hence,

$$\begin{aligned} p_1^*=125,\quad p_2^*=175, \end{aligned}$$

with associated utilities being:

$$\begin{aligned} {\hat{u}}_1(x^*)=u_1(p^*)=875,\quad {\hat{u}}_2(x^*)=u_2(p^*)=325, \end{aligned}$$

and the migration costs on the routes: \(C_{r_1}=0\), \(C_{r_2}=0\), \(C_{r_3}=300\), \(C_{r_4}=162.50\). Moreover, \(\lambda _1^*=\lambda _2^*=0\). The optimal Lagrange multiplier associated with the regulation constraint is: \(\mu _1^*=387.50\). One can see that route \(r_3\) is too expensive and will not be used under the refugee migratory flow pattern.

Now, refugees in country 1, at equilibrium, enjoy a higher utility than before the regulation was imposed of 875 versus 775, an increase of about 13%. On the other hand, refugees in country 2 now experience a lower utility, than before the regulation was imposed. They are no longer all free to migrate because of the imposed regulatory upper bound of 25 limiting the migration from country 2 to country 1. Under the regulation, those who manage to migrate enjoy a higher utility than before of 425, but those who are left behind in country 2 experience a lower utility than in the case without regulations—of 325, revealing a drop of about 30%.

4 Supernetwork transformation to a traffic network equilibrium problem with fixed demands

Before presenting the computational procedure that we will apply to compute solutions to larger refugee migration problems, we establish a supernetwork equivalence between the refugee model with multiple paths, but with a single class, and a traffic network equilibrium problem with fixed demands. Such a connection has never been made before and enables the transfer of a plethora of algorithms from transportation science to solve refugee migration problems. We note that the multiclass transformation into a single class/mode problem readily follows by applying standard multicopy techniques as proposed by Dafermos [6].

Specifically, recall, cf. Fig. 1, that, in the refugee network, there are origin nodes where the initial populations, including the refugees, are located. Moreover, in order to capture that some of the refugees may remain at their origin locations, we also have as many destination nodes in the original network as there are origin nodes. Plus, in order to reflect that a volume of the population may remain at each location, we have a path consisting of a single link joining an origin node to its “companion” destination node, with an associated cost, as noted earlier, of zero. Also, in the original network there are paths/routes consisting of one or more links joining origin nodes to destination nodes reflecting possible routes that the refugees may take; see, for example, multiple illustrations of this in the real world on different continents on actual maps in Conant [5].

We now construct the supernetwork representation—cf. Fig. 3. We add a single “super” destination node D to the bottom of the original network in Fig. 1 and a link joining each original destination node to node D. The origin/destination (O/D) pairs on the supernetwork, which will be equivalent to a traffic network equilibrium problem with fixed demands (cf. Dafermos [7], Smith [48]) are now as follows. There are as many O/D pairs in the supernetwork as there are origin nodes, that is, n, with the O/D pairs consisting of the following pairs of nodes: \(w_1=(1,D)\), \(\ldots \), \(w_n=(n,D)\). Note that each origin node in the supernetwork corresponds to the respective origin node in the original network. The fixed demands associated with these O/D pairs, denoted by \(d_{w_1}\), \(d_{w_2}\), \(\ldots \), \(d_{w_n}\) are set, respectively, equal to: \({\bar{p}}_1\), \({\bar{p}}_2\), \(\ldots \), \({\bar{p}}_n\). The costs on the links on the supernetwork that are the same as in the original network remain the same. The costs on the new links joining destination nodes 1, \(\ldots \), n, with node D capture the utility information at the respective destination and are assigned, respectively, as: \(-u_1(p)\), \(\ldots \), \(-u_n(p)\).

Fig. 3
figure 3

Supernetwork transformation of the sample refugee network topology

It is straightforward to verify that the flows on paths connecting each O/D pair on the supernetwork will coincide with the respective path flows on the original network and the conservation of flow equations (2) and (4) will hold. Furthermore, the final populations (as in (3)) at the destination locations coincide with the link flows on the respective additional links in the supernetwork in Fig. 3, terminating in node D.

We now prove that the equilibrium conditions according to Definition 1 are also satisfied by the TNE flow on the supernetwork in Fig. 3 with the additional link costs as defined above.

Observe that, from the construction of the supernetwork and the associated link costs, in equilibrium, if the flow on the path consisting of the original direct link from an origin node i to its companion node and the super demand node D is positive, then the cost on that path is precisely equal to minus the utility, \(-u_i(p^*)\). Moreover, if the flow on any other path from that same origin node to D is positive, and, say, that it goes through a node j to node D, then the cost on that path is equal to \(C_r-u_j(p^*)\), where \(C_r\) is the cost on the original such path r and \(-u_j(p^*)\) is the cost on the final link to node D. Because of the structure of the supernetwork, we know that these two costs are equalized, in the traffic network equilibrium, and, hence, we have that:

$$\begin{aligned} C_r-u_j(p^*)=-u_i(p^*), \end{aligned}$$

or, equivalently:

$$\begin{aligned} u_i(p^*)+C_r=u_j(p^*), \end{aligned}$$

which is precisely the first part of the equilibrium condition (6) (since this holds true for any such paths and nodes). Clearly, the second part of (6) also holds since any more expensive path connecting each O/D pair will not be used. The second part of the equilibrium conditions, (7), is easy to see will hold, as well, if the initial population at an origin location node is insufficient, in which case the corresponding \(\lambda _i^*\) may be positive.

Also, it is well-known that the path flow variational inequality formulation of the fixed demand traffic network equilibrium problem (cf. Dafermos [7] and Nagurney [26]) is of the form: determine \({\hat{X}}^*\in {{{\mathcal {K}}}}^1\) such that

$$\begin{aligned} \langle {\hat{C}}({\hat{X}}^*),{\hat{X}}-{\hat{X}}^*\rangle \ge 0,\quad \forall {\hat{X}}\in {{{\mathcal {K}}}}^1, \end{aligned}$$
(17)

where \({\hat{C}}\) is the vector of path costs, \({\hat{X}}\) is the vector of path flows, and \({{{\mathcal {K}}}}^1\) is the feasible set consisting of the nonnegativity assumption on the path flows and with the demands for the O/D pairs being satisfied by the path flows. With the path flows and the path costs as defined above on the supernetwork, along with the feasible set, it is clear that variational inequality (17) coincides with variational inequality (15).

The supernetwork reformulation is relevant even for the refugee model with regulations since the regulations become side constraints for the former and variational inequality (17) still applies but with the feasible set \({{{\mathcal {K}}}}^2\) including also the side constraints. We assume that \({{{\mathcal {K}}}}^2\) is nonempty.

We provide the full supernetwork equivalence for the numerical example(s) in Sect. 6.

Through the supernetwork equivalence, one can avail oneself of many algorithms that have been applied to solve traffic network equilibrium problems (cf. Bar-Gera [1], Nagurney [22, 23], Patriksson [43], Sheffi [47], and the references therein), which have, also, because of the longer history of that area of application, been used effectively and efficiently to solve large-scale real world transportation networks in cities consisting of thousands of nodes, links, and origin/destination pairs of nodes (see Perederieieva et al. [44] and the references therein).

5 The algorithm and numerical examples

Our numerical examples are inspired by the refugee flows from Mexico to the United States, an issue that has been receiving a lot of attention in the press. Specifically, the baseline network for the numerical examples is depicted in Fig. 4. The top nodes correspond to origin nodes, representing cities/towns, with the first four nodes being locations in the southwestern United States and the next three top nodes corresponding to locations in Mexico, where refugee origin nodes have been identified. The bottom nodes in the network in Fig. 4 are the destination node locations. We consider a single class of refugee.

Fig. 4
figure 4

The baseline refugee network topology for the numerical examples

We consider four examples: Example 1 through Example 4, with Example 1 being the example without regulations.

Example 1

The initial populations at the origin nodes are:

$$\begin{aligned} {\bar{p}}_1= & {} 1,400,000,\quad {\bar{p}}_2=20,000,\quad {\bar{p}}_3=70,000, \\ {\bar{p}}_4= & {} 260,000,\quad {\bar{p}}_5=50,000,\quad {\bar{p}}_6=225,000,\quad {\bar{p}}_7=30,000. \end{aligned}$$

The utility functions associated with these locations are as follows:

$$\begin{aligned} u_1(p)= & {} -p_1+3,000,000,\quad u_2(p)=-2p_2+200,000,\quad u_3(p)=-p_3+1,500,000, \\ u_4(p)= & {} 3p_4+900,000,\quad u_5(p)=-p_5+100,000,\quad u_6(p)=-p_6+300,000, \\ u_7(p)= & {} -p_7+100,000. \end{aligned}$$

From the above utility functions, one can see that the locations in the United States are more attractive than those in Mexico, due to the significantly larger fixed utility term in the corresponding utility functions.

We now provide the link cost data. Please refer to the network in Fig. 4 for the link locations.

Specifically, the link costs associated with remaining at one’s location at nodes 1 through 7, respectively, are, as discussed earlier, all equal to 0.00, and, hence, we have that:

$$\begin{aligned} c_1(f)=c_2(f)=c_3(f)=c_4(f)=c_{11}(f)=c_{15}(f)=c_{18}(f)=0.00. \end{aligned}$$

The costs associated with the refugee migrations are, in turn, as follows:

$$\begin{aligned} c_5(f)= & {} .00006f_5^4+6f_5+4f_6+200,\quad c_6(f)=7f_6+3f_8+300,\\ c_7(f)= & {} .00008f_7^4+8f_7+2f_8+400, \\ c_8(f)= & {} .00004f_8^4+5f_8+2f_{10}+450,\quad c_9(f)=.00001f_9^4+6f_9+2f_{10}+300,\\ c_{10}(f)= & {} 4f_{10}+f_{12}+400, \\ c_{12}(f)= & {} 8f_{12}+2f_{13}+100,\quad c_{13}(f)=.00001f_{13}^4+7f_{13}+3f_9+50,\\ c_{14}(f)= & {} 8f_{14}+3f_9+100, \\ c_{16}(f)= & {} 3f_{16}+f_{12}+100,\quad c_{17}(f)=.00003f_{17}+3f_{17}+50. \end{aligned}$$

Observe that the above cost functions are nonseparable and several have flow links raised to the fourth power. The former characteristic allows one to model interactions whereas the latter captures increased congestion and, hence, cost.

One can see from the network in Fig. 4, that the refugees located at location node 5 in Mexico have three available paths to reach the United States, whereas those at node 6 and node 7 have a single path each, all with associated costs on the routes, which consist of multiple links. Those residing at location nodes 1 through 4 in the United States are not interested in migrating to Mexico. Of course, those living at nodes 5, 6, and 7 can also elect to remain in Mexico.

The routes are defined as consisting of the following respective links and enumerated as follows:

$$\begin{aligned} r_1= & {} (1),\quad r_2=(2),\quad r_3=(3),\quad r_4=(4), \\ r_5= & {} (5,6),\quad r_6=(7,8),\quad r_7=(9,10),\quad r_{8}=(11), \\ r_9= & {} (12,13,14,10),\quad r_{10}=(15), \\ r_{11}= & {} (16,17),\quad r_{12}=(18). \end{aligned}$$

For the solution of the problem, as noted earlier, we first construct the supernetwork equivalence with the supernetwork topology as in Fig. 5 and the O/D pairs, the links and paths, the link costs on the new links, the demands, and the path costs, defined accordingly.

Fig. 5
figure 5

The supernetwork topology for the baseline example

We implemented the Euler method, embedded with the exact equilibrium algorithm as described in Nagurney and Zhang [40]; see also the book by Nagurney and Zhang [39]. We initialized the algorithm as follows. The initial populations were equally distributed among all the paths. The sequence \(\{a_\tau \}\) in the Euler method satisfied the conditions required for convergence and was set to: \(.1\{1,{{1}\over {2}},{{1}\over {2}},{{1}\over {3}}, {{1}\over {3}}, {{1}\over {3}}, \ldots \}\). The algorithm was deemed to have converged if the absolute value of the difference between each successively computed path flow differed by no more than \({10}^{-7}\). The computer system utilized was a Linux system at the University of Massachusetts Amherst.

The computed equilibrium path flow pattern for Example 1 is:

$$\begin{aligned} x_{r_1}^*= & {} 1,400,000.00,\quad x_{r_2}^*=20,000.00,\quad x_{r_3}^*=700,000.00,\quad x_{r_4}^*=260,000.00, \\ x_{r_5}^*= & {} 400.25,\quad x_{r_6}^*=336.64,\quad x_{r_7}^*=317.17,\quad x_{r_8}^*=48,945.95, \\ x_{r_9}^*= & {} 290.25,\quad x_{r_{10}}^*=224,709.75,\quad x_{r_{11}}^*=199.54,\quad x_{r_{12}}^*=29,800.46, \end{aligned}$$

which corresponds to the following equilibrium link flow pattern:

$$\begin{aligned} f_1^*= & {} 1,400,000.00,\quad f_2^*=20,000.00,\quad f_3^*=700,000.00,\quad f_4^*=260,000.00, \\ f_5^*= & {} 400.25,\quad f_6^*=400.25,\quad f_7^*=336.64,\quad f_8^*=336.64, \\ f_9^*= & {} 317.17,\quad f_{10}^*=607.41,\quad f_{11}^*=48,945.95,\quad f_{12}^*=290.25, \\ f_{13}^*= & {} 290.25,\quad f_{14}^*=290.25, \quad f_{15}^*=224,709.75,\quad f_{16}^*=199.54, \\ f_{17}^*= & {} 199.54,\quad f_{18}^*=29,800.46. \end{aligned}$$

The computed equilibrium populations at the four locations in the United States and at the three locations in Mexico are:

$$\begin{aligned} p_1^*= & {} 1,400,736.88,\quad p_2^*=20,607.42,\quad p_3^*=700,000.00, \\ p_4^*= & {} 260,199.55,\quad p_5^*=48,945.95,\quad p_6^*=224,709.75,\quad p_7^*=29,800.46. \end{aligned}$$

Note that these equilibrium location populations correspond, respectively, to the equilibrium link flows on links 19 through 25 in the supernetwork in Fig. 5.

The associated incurred utilities at the equilibrium at the locations are:

$$\begin{aligned} u_1(p^*)= & {} 1,599,263.13,\quad u_2(p^*)=158,785.17,\quad u_3(p^*)=800,000.00, \\ u_4(p^*)= & {} 119,401.38,\quad u_5(p^*)=51,054.05,\quad u_6(p^*)=75,290.25,\quad u_7(p^*)=70,199.55. \end{aligned}$$

We also, for completeness, report the path costs at the computed equilibrium flows since it is easy to then verify that the equilibrium conditions, as stated in Definition 1, hold. Specifically, the incurred path costs at the equilibrium are:

$$\begin{aligned}&C_{r_1}(x^*)=C_{r_2}(x^*)=C_{r_3}(x^*)=C_{r_4}(x^*)=0.00, \\&\quad C_{r_5}(x^*)=1,548,207.50,\quad C_{r_6}(x^*)=1,548,196.38,\quad C_{r_7}(x^*)=107,729.39, \\&\quad C_{r_8}(x^*)=0.00,\quad C_{r_9}(x^*)=83,500.69,\quad C_{r_{10}}(x^*)=0.00, \\&\quad C_{r_{11}}(x^*)=49,200.66,\quad C_{r_{12}}(x^*)=0.00. \end{aligned}$$

We now proceed to impose regulations and construct Examples 2 through 4, which are analogues of Example 1, but with regulations.

Examples 2, 3, and 4

In Example 2, we considered the following scenario: The US government has told the Mexican government that it is restricting the flow on route \(r_5=(5,6)\) to zero; essentially resulting in the elimination of this path, since its processing facilities are experiencing delays due to congestion. The rest of the data remain as in Example 1.

In Example 3, we studied the following scenario: Route \(r_6=(7,8)\) is now unavailable to refugees, but the other routes and data remain as in Example 1.

And, in Example 4, we investigate the scenario that the United States is concerned about the influx of refugees and both routes \(r_5\) and \(r_6\) from Mexico are, in effect, banned/eliminated.

In order to enable cross comparisons with the examples with regulations and with the baseline Example 1, in Table 2, we report the computed equilibrium path flows and in Table 3 we report the associated path costs for all four examples. In Table 4, we report the associated computed equilibrium populations, whereas the incurred utilities at the locations are reported in Table 5.

Table 2 Equilibrium path flows for Examples 1–4
Table 3 Equilibrium path costs for Example 1–4
Table 4 Equilibrium populations at the locations for Example 1–4
Table 5 Equilibrium utilities at the locations for Example 1–4

One can see, from Table 5, that (as occurred also in the illustrative examples in Sect. 3), under the regulations (as in Examples 2 through 4), the utility of those subject to regulations, as in the case of those in node 5, is reduced. On the other hand, those in location 1, which now has a lower flow of refugees, experience a higher utility. Also, with both refugee routes blocked to the United States, the population that remains at node 5 is the highest in Example 4, as compared to the value in Example 2 and Example 3.

We remark that although the examples are stylized, they nevertheless, contain features representative of refugee networks in reality.

6 Summary, conclusions, and suggestions for future research

In this paper, we constructed refugee migration network equilibrium models that incorporate multiple classes of refugees, along with routes from origin to destination nodes that can consist of multiple links that capture congestion. Such features have been noted in practice but, heretofore, not captured in models of human migration. Moreover, we utilized a user-optimization perspective, rather than a system-optimization perspective. In the former, the refugees/migrants seek locations to migrate to that are optimal from an individual perspective as opposed to resettling refugees in a system-optimized manner.

In addition, we established that the more general model, for which we also proposed constraints associated with governmental migration regulations, could be transformed, using a supernetwork construct, into a traffic network equilibrium with fixed demands. This connection expands the theory of network equilibrium problems and allows for the transferral of algorithms that have been developed over many years for the solution of traffic network equilibrium problems to problems of refugee and human migration. In addition, we noted that regulations in this context can be viewed as side constraints associated with traffic networks.

We provided both small numerical examples, in order to illustrate the modeling framework, as well as larger scale examples, inspired by the refugee crisis from Mexico to the United States. For the solution of the latter numerical examples, we utilized the Euler method, embedded with the exact equilibrium algorithm. Our examples show who may benefit from the imposition of certain regulations and who may lose. We expect that research in the future may investigate additional types of regulations of refugee networks, for which other computational schemes may be utilized, with further exploitation of network structure.