A new perspective for community detection

Community detection in a network plays a crucial role in the economic and financial contexts, specifically when applied to the World Trade Network. We provide a new perspective in which clusters of strongly interacting countries are identified by means of a specific distance criterion. We refer to the Estrada communicability distance and the vibrational communicability distance, which turn out to be particularly suitable for catching the inner structure of the economic network. The methodology is based on a varying distance threshold and it is effective from a computational point of view. It also allows an inspection of the intercluster and intracluster properties of the resulting communities. The numerical analyses highlight peculiar relationships between countries and provide a rich set of information that can hardly be achieved within alternative clustering approaches.

Introduction
International trade is based on a set of complex relationships between different countries. Both connections between countries and bilateral trade flows can be modelled as a dense network of interrelated and interconnected agents. A longstanding problem in this field is the detection of communities, namely subset of nodes among which the interactions are stronger than average. Indeed, the community structure of a network reveals how it is internally organized, highlighting the presence of special relationships between nodes, that may not be caught from direct empirical analyses. In other words, the community structure refers to the occurrence of clusters of nodes that are more densely connected internally than with the rest of the network. In this framework, a specific role is assumed by the distance between nodes. For instance, the neighbours of a given node are immediately connected to such a node and they can affect its status in the most direct way. Nonetheless, more distant nodes can influence this node while passing through intermediary ones. In the economic field, a network perspective is indeed based on the idea that indirect trade relationships may be important (see, e.g., [31]). The authors in [1] explain the impact of shocks on a given country by indirect trade links. Based on a global VaR approach, [16] shows that countries that do not trade (very much) with the U.S. are largely influenced by its dominance over other trade partners linked with the U.S. via indirect spillovers. In [62], the bilateral trade is assumed not independent of the production, consumption, and trading decisions made by firms and consumers in third countries. A measure of the distance between nodes that also considers indirect connections is therefore crucial to catch deep interconnections between nodes. In this work, we will focus on two measures of distance or metrics on the network: the Estrada communicability distance [22] and the vibrational communicability distance [23]. They both go beyond the limits of the immediate interaction between neighbours and they look simultaneously, albeit differently, at all the possible channels of interactions between nodes. The nearest two nodes are in each metric, the stronger is their interaction or, in other words, the higher is the level of communicability between them. With this paper we contribute to the literature proposing a new method that exploits such metrics to inspect the mesoscale structure of the network, in search for strongly interacting groups of nodes. Indeed, our purpose is to provide a new methodology for detecting the community structure of a network exploiting these suitable distances. In other words, we aim at finding cluster of nodes by means of a metric able to catch both direct and indirect relationship between countries. In particular, as in the classical procedure provided by Girvan and Newman [46] we search for the maximum modularity, that, as well-known, is a way to measure if a specific mesoscopic description of the network in terms of communities is more or less accurate. Differently from [46], we group together nodes whose interaction is lower than a specific threshold. The best partition will be identified by the threshold that maximizes the modularity.
Our proposal provides several advantages representing a viable alternative to classical methodologies for community detection. Firstly, the method is very efficient from a computational viewpoint. Indeed, given the specific distance matrix, the optimal solution can be easily evaluated varying the threshold. Classical Girvan and Newman methodology is instead a NP-hard problem due to the fact that the space of possible partitions grows faster than any power of the system size. For this reason, several heuristic search strategies have been provided in the literature to restrict the search space while preserving the optimization goal (see, e.g., [2], [17], [10], [45]). Secondly, we cluster nodes going beyond the interactions between neighbours and considering all possible channels of interaction between nodes. Thirdly, we allow for a degree of flexibility by introducing a threshold. Varying the threshold, it is possible to depart from the optimal solution so that only the strongest (or the weakest) channels of communications emerge. Finally, the procedure offers a set of indicators that allow to exploit main characteristics of the communities detected as well as the relevance of countries inside the community and in the whole network. The theoretical model is validated through empirical experiments based on the World Trade Network (WTN). Results confirm a number of intuitions in mainstream economic literature. It is confirmed the relevance of physical variables such as geographical distance (see, e.g., [4]). Geographical proximity matters and trades are to some extent regionalized. These results are also in line with [47] that shows that existing preferential trade agreements are not strongly distorting the geography of trade patterns, at least at the aggregate level. Two main communities are detected. Top community is characterized by North-America, Asian countries and Australia. In line with [65], emphasizing tight links, we obtain that although the strong trade relationships with USA and Germany, China became regionally attractive and restored the leadership of Asia-Oceania community. European community is highly centralized around founding members of the European Economic Community with the central role of Germany. High income countries in Northern Europe are instead in a separate community with a less relevant role in the network. The paper is organised as follows. After a short review of the literature in Section 2, main preliminaries and the definitions of the communicability functions are revised in Section 3. These functions lead to two important metrics on networks, which are described in Section 4. Section 5 contains the description of the proposed methodology, which is also tested on a suitable toy-model. In Section 6, we test our methodology to the World Trade Network. In particular, main characteristics of the network are described in Section 6.1. Main steps of the methodology are summarized in Section 6.2. We report in Sections 6.3 and Subsection 6.3.1 main results based on communicability and resistance distance, respectively. We show how the proposed methodology is able in capturing key economic clusters as well as in providing additional insights in terms of intercluster characteristics and of countries' relevance in the community and in the whole network. Conclusions follow. Technical details are left in Appendices A and B.

Literature Review
Community detection is an important topic in the analysis of the topological structure of complex systems. This is confirmed by the fact that its importance has grown over time in light of the remarkable progress in the description of large networks, together with the development of new powerful data analysis tools. These advances have made it possible to extend the field of applicability of the theory not only to networks of enormous dimensions but also to weighted networks and direct networks [9,12,27,49]. The definition of what a community is within a network is not univocal but there is a broad general consensus in believing that it maximizes a modularity function according to Newman's definition [46].
More recently the role of non-local interactions between nodes has been highlighted, that is interactions that do not exclusively involve the immediate neighbours of a given node. In particular, results connected to the idea of communicability introduced by Estrada in 2004 have proved to be extremely effective [19,21,22,25]. All the more so by allowing a metric different from the shortest path metric to be introduced on the network. The purpose of this new metric is precisely to take into consideration long-range interactions between institutions. Some important similarities can be found between this new metric and the resistance distance, a well-known metric in network theory derived from the study of electric circuits [23,40,43], and its interpretation in terms of vibrational communicability [3,6,24,59].
An area in which these concepts allow us to gain a deep insight into the hidden structures of the network is properly the WTN. The topology of the world trade web has been extensively analysed over time [28, 33-35, 44, 53]. The behaviour of international trade flows, the impact of globalization on the international exchanges, the presence of a core-periphery structure or the evolution of the community centres of trade, are just some of the issues addressed by the recent developments [5,14,29,54,57]. Many works have dealt with the network from a multi layers perspective [4,56] or aim to emphasize financial implications of the world trade or contagion processes on the network [7,8,15,30,32,36,50,52,60,63].
For instance, centrality of a given institution as best spreader node in a contagion process has been discussed by [43,59], who define the avalanche size of a node as the number of subsequently collapsed institutions starting from a given institution's collapse. In that way they are able to identify the most dangerous crisis epicentre of the network, be it a country or a financial institution.
The impact of topology and metric properties on the stability and resilience of an economic or financial system has been widely studied in order to describe the large-scale pattern of dynamical processes inside the network [39,48,55]. These processes determine the subsequent diversification of the export of a country, which can be compared with descriptive empirical indices of its potential growth, such as the one introduced in a very fruitful way in [38].

Communicability in complex networks
The idea of communicability on a network is based on the ways in which a pair of nodes can communicate, namely through walks connecting them. In the literature, two different definitions of communicability have been introduced: the Estrada Communicability and the Vibrational Communicability ( [21,24]). We recall them in this section.
First of all, we briefly remind some preliminary definitions. A network is formally represented by a graph G = (V, E) where V and E are the sets of n nodes and m edges, respectively. Two nodes i and j are adjacent if there is an edge (i, j) ∈ E connecting them. The network is undirected if both (i, j) and (j, i) are elements of E. A i − j-path is a sequence of distinct vertices and edges between i and j. The shortest path, or geodesic, between i and j is a path with the minimum number of edges. The length of a geodesic is called geodesic distance or shortest path distance Adjacency relationships are represented by a binary symmetric matrix A (adjacency matrix). Graphs considered here will be always connected and without loops; in this case a ii = 0 ∀i = 1, ..., n. We denote with λ 1 ≥ λ 2 ≥ · · · ≥ λ n the eigenvalues of A, and ϕ i , i = 1, ..., n the corresponding eigenvectors. The degree k i of a node i is the number of edges incident on it. The diagonal matrix whose diagonal entries are k i is K. The Laplacian matrix is L = K − A. L is positive semidefinite symmetric. We denote the eigenvalues of L by µ 1 ≥ µ 2 ≥ · · · > µ n = 0 and ψ i , i = 1, ..., n the corresponding eigenvectors.
A graph G is weighted when a positive real number w ij > 0 is associated with the edge (i, j). We define the strength s i as the sum of the weights of the edges adjacent to i. The definition of geodesic path still holds, and it is a weighted path with the minimum sum of edge weights. In this case, the adjacency matrix is a non-negative symmetric matrix W. When w ij = 1 if (i, j) ∈ E, then the graph is unweighted. Thus, the unweighted case can be viewed as a particular weighted one. We will use in the paper the general concept of weighted graphs and we denote a weighted graph with its weights simply as weighted network.
A weighted network is directed if every edge is directed, i.e. if (i, j) is an element of E, not necessarily (j, i) is in E. The weighted adjacency matrix W is not symmetric. The in-strength s in i and the out-strength s out i of a node i are, respectively, the sum of the strengths of all the arcs pointing to or starting from such a node: where W i represents the i-th row of the matrix W and u ∈ R n the vector whose components are all 1.

Estrada Communicability
The Estrada communicability [21] between two nodes i and j is defined as: As the ij-entry of the k-power of the adjacency matrix A counts the number of walks of length k starting at i and ending at j, G ij accounts for all channels of communication between two nodes, giving more weight to the shortest routes connecting them. It can also be interpreted as the probability that a particle starting at i ends up at j after wandering randomly on the complex network. The communicability matrix is denoted by G.
By definition, it follows that G ij > 0. Moreover, G ij can be conveniently expressed using the spectral decomposition of A as follows [21]: where ϕ k (i) is the i-component of the k-th eigenvector associated with λ k . It is worth noting that since G ii characterizes the importance of a node according to its participation in all closed walks starting and ending at it, we recover the so-called subgraph centrality (see [25]).
In the case of an undirected weighted network the communicability function is defined as where S is the diagonal matrix whose diagonal entries are the strengths of the nodes. We will call this weighted communicability. The weighted communicability is particularly suitable to be applied to the study of input-output networks.
The definition of communicability can be generalized in order to describe situations in which an external factor may affect the structure of connections. Indeed, complex networks are continuously exposed to external stresses, which are independent of the network architecture, but that can affect its structure and its evolution. To take into account the effect of the external environment on the network communicability, a real parameter β can be introduced as a weight-gauging factor which multiplies the original weights assigned to each link, i.e. by multiplying the adjacency matrices A or W by the constant β. See [22] for a more detailed explanation about this generalization.

Vibrational Communicability
Communicability can be alternatively defined through the following model from Physics. Let us suppose that nodes of the network are objects of negligible identical mass connected by springs in a plane grid. Nodes can oscillate in the direction perpendicular to the plane and the displacement of the node i from its rest position is z i . The elastic force applied to node i is given by where K is the common elastic constant of each spring. An elastic potential energy can be assigned to each perturbed spring and the potential energy of all the springs connected with node i is given by The overall potential energy of the network is therefore 1 where L ij is the ij-entry of L.
The reciprocal influence of two nodes i and j in their positions z i and z j is computed by means of the Green's function, according to the classical Boltzmann's distribution [23,24]. This mutual influence can be interpreted as the correlation function between the displacements z of two nodes in the network: where β is a constant and Z = e −βU dz is the partition function. Using the non-zero eigenvalues of L, Z can be expressed as (see Appendix 2) so that the correlation function can be rewritten in the final form (see Appendix 2): where ψ k is the eigenvector associated with µ k . Introducing the Moore-Penrose pseudo-inverse of the Laplacian L + [6,37], the vibrational communicability between nodes i and j is defined as The vibrational communicability matrix is denoted by G v . In the remainder of the paper we will assume β = 1 and K = 1, so that G v ij = L + ij .

Metrics on networks
Metric properties play an important role in the study of the structure and dynamics of networks. The best known metric is the so-called shortest path distance. In the literature other metrics have been defined, each one stressing different features of the network. We remind the definitions of communicability distance and resistance distance, in view of their following application to the WTN.

Communicability Distance
The communicability distance ξ ij is defined as (see [19]): As already observed, G ii is the subgraph centrality of i and it measures the amount of information that starts from and returns to node i after having wandered through the network. On the other hand, G ij measures the amount of information transmitted from i to j. Notice that the word information is meant in its broadest sense. Therefore, information flow can be any kind of flow along edges: money, current, traffic and so on. Thus, the quantity ξ ij accounts for the difference in the amount of information that returns to the nodes i and j and the amount of information exchanged between them. The greater is G ij , the larger the information exchanged and the nearer are the nodes; the greater are G ii or G jj , the larger the information that comes back to the nodes and the farther are the nodes. In a matrix form, ξ ij can be expressed as follows: , . . . , G nn ] T is the vector of subgraph centralities. Since ξ ij is a metric, then G ii + G jj ≥ 2G ij , i.e., no matter what the structure of the network is, the amount of information absorbed by a pair of nodes is always larger than the amount of information transmitted between them.

Resistance Distance
The vibrational communicability distance between i and j is defined as (see [23,59]): Formula 10 can be written in a more suitable way. Indeed, recalling that G v ij = L + ij , we have: where e k , k = 1, . . . , n, is the standard basis in R n and J = uu T is the matrix whose entries are all 1. Note that in the previous chain of equalities we made use of the following expression of the pseudo-inverse [37]. Equation 11 offers an interesting interpretation in the networked system. We synthesize here the idea of this interpretation, referring to the Appendix B for a more detailed discussion. Let v = [v 1 , v 2 , . . . , v n ] T be a vector representing the attributes of the nodes -for instance, the Gross Domestic Product (GDP) of a Country or the assets of a financial institution -and suppose the existence of a flow of currents in the network. The operator L + 1 n J −1 allows to obtain the state vector on the nodes starting from this flow vector. In our case, we want a flow equal to +1 from node i, a flow equal to −1 into node j and a flow equal to 0 for the other ones. The inner product with (e i − e j ) of formula 11 gives v i − v j , namely, the difference between attributes of nodes i and j. This gradient produces exactly the flow +1 from node i and −1 to node j. If v i − v j is big, we need a big difference in order to produce such a unit flow and so we have a big resistance between nodes i and j. If v i − v j is small, it is enough a low difference in order to produce such a unit flow and so we have a low resistance between nodes i and j. If ω ij is big we have a high resistance distance between i and j. Therefore, these two nodes do not communicate easily. Vice versa a low value of ω ij means a high level of communication between the nodes. ω ij is called effective resistance between nodes i and j. and Ω = [ω ij ] is the resistance matrix.
In the literature, it is known an important close form for L + in terms of Ω: which allows us to rewrite the diagonal elements of the matrix L + in a useful form 2 is the effective graph resistance (or Kirchhoff index ) of the network, i.e. the sum of the resistances between all possible pairs of nodes in the graph (see, e.g., [40]). R reflects the overall transport capability of the network: the lower R, the better the network conducts flows. In particular, it has been shown that this index is able to catch the vulnerability of a connection between a pair of nodes and, therefore, it is a suitable tool for assessing the ability of a network to well react when it is subject to failure and/or attack (see [11,18,61]).
Additionally, the relevance of a node can be quantified by exploiting the effective resistances. Indeed, the best spreader (or best connected) node in the network is the node i that minimizes the the quantity n j=1 ω i j = (Ωu) i , i.e. the sum of all resistances between itself and any other node in the network. Since L + ii equals the difference between the average resistance between node i and all the other nodes in the network and the overall network's mean resistance, then the best spreader node i is the one such that L + i i ≤ L + jj for any j = i . Node i can be regarded as the best diffuser of a flow to the rest of the network, and, to some extent, it is the most influential with respect to a diffusion process inside the network, since it guarantees the highest flow toward other nodes (see [59]). The reciprocal of L + ii can then be regarded as a centrality measure of a node and it is called vibrational centrality. Best diffuser means that most of the information coming out from this node is absorbed by other nodes. If L + ii is big, then most of this information comes back to node i and doesn't reach other nodes.

The model
As discussed in the previous section, represent the two metrics induced on the network by the Estrada communicability and the vibrational communicability, respectively.
In an economical context, referring to the international trade network, they measure how well two countries, or companies, communicate in terms of commercial and trade exchanges. For instance, the attributes on nodes may be identified with the GDP and the currents along nodes with the total trade or money flow between two countries. Information on the network may be replaced by money flow. Therefore the quantity ξ ij of equation 9 accounts for the difference in the amount of money flow that returns to the nodes i and j and the amount of money flow exchanged between them. The bigger is G ij , i.e. the money flow exchanged, the nearer are the nodes; the bigger are G ii or G jj , i.e. the amount of money flow that comes back to the each node, the farther they are.
A similar interpretation holds for ω ij . In a trade network ω ij accounts for the difference between the mean resistance to export a given money flow from each country and the correlation between them. The bigger is G v ij , the more interconnected they are and the nearer they are in the resistance metric; the bigger are G v ii and G v jj , the more isolated they are in the network and between them and the farther they are.
We formulate our proposal in light of these observations. In what follows, we will refer to the communicability distance ξ, but similar arguments may be repeated identically for the resistance distance ω.
We consider members of the same cluster nodes whose mutual distance is below a given threshold ξ 0 . Specifically, we construct a new community graph where the elements of the adjacency matrix M = [m ij ] are given by: , being ξ min and ξ max the minimum and the maximum distances between couples of nodes, respectively. In this way, clustered groups of nodes that strongly communicate emerge, in dependence of the threshold. If ξ 0 is high enough, all nodes in the network are at a mutual distance lower than the threshold and the whole network behaves like a unique community. As ξ 0 decreases, there will be nodes too far, such that to be considered disconnected or members of different clusters, and we witness the emergence of islands of still connected nodes.
Hence, the number of communities depends on the threshold, precisely it increases as ξ 0 decreases. It is important to observe that the proposed method cannot choose any a priori optimal number of communities. Rather, as the classic Girvan-Newman approach [46], the optimal partition, and the corresponding threshold, are chosen according to the following optimization problem for the modularity Q as a function of ξ 0 : where x ij are binary variables equal to 1 if two nodes are in the same cluster and 0 otherwise, and S is total strength of the network. It is worth noting that the modularity function is computed on a nodes partition emerging through the community graph, in dependence on the threshold ξ 0 . The best partition of the network in communities will be identified with the one which maximises the modularity function Q. The pseudocode defining the proposed methodology will be better described in the next section.

An illustrative example
We start by testing our methodology on a simple example. Let us consider the weighted undirected network in Figure 1. The network has 10 nodes and 32 edges. The thickness of links is proportional to weights. The network allows to easily identify two natural communities, which are highlighted by the two closed lines containing nodes 1 to 5 (on the left) and nodes 6 to 10 (on the right). We compute the Estrada communicability matrix G, then we get the communicability distance matrix Ξ. The nearest nodes are 1 and 3 with a communicability distance equal to ξ min = ξ 13 = 1.18 and farthest nodes are 3 and 6 with a communicability distance equal to ξ max = ξ 36 = 1.49. Figure 2 summarizes the number of communities identified for alternative thresholds. The blue line represents the number of communities while the red line represents the value of modularity of the corresponding partition. When the threshold is equal to 1.38 all nodes are connected and the network is partitioned in a single community with modularity Q = 0. As the threshold decreases below 1.38, the network begins to split into disconnected components. When the threshold becomes lower than the minimum distance, the network is partitioned into ten communities, namely, each node belongs to a different community. The best partition according to the maximum modularity criterion splits the network into two clusters, which are easily identified with the two expected natural communities. The composition of the communities for alternative thresholds is reported in Figure 3. It is noticeable that, lowering the threshold, the procedure allows to disentangle tightest relationships. For instance, when ξ 0 = 1.23 only nodes connected by edges with highest weights are kept in the same community. Similar results are derived applying the procedure based on the vibrational communicability. The nearest nodes are 2 and 5 with a resistance distance equal to ω min = ω 25 = 0.29 and farthest nodes are 1 and 6 with a resistance distance equal to ω max = ω 16 = 0.42. Again if we move the threshold from the maximum distance to the minimum distance, we get an increasing number of communities from 1, the whole network, to 10, isolated nodes. The best partition according to the maximum modularity criterion splits the network into the two expected communities, as shown in Figure 4.

Application to the World Trade Network
In this Section, we apply the proposed model in order to detect relevant communities of countries in the WTN. As described before, the method aims at grouping countries on the basis of their distance in the network. Two alternative distance functions will be tested. On the one hand, we find clusters exploiting communicability distance. Therefore we detect how much two countries are close in the network considering all possible weighted walks between them. On the other hand, we select clusters by means of resistance distance. In this case countries are grouped together if they have a similar relevance in the network in terms of vibrational centralities as well as if they are correlated in terms of their expositions towards common countries. We start with a general description of the dataset and we provide main characteristics of the WTN. Then, we briefly summarize the main steps of the methodology also providing a pseudo-code. Finally, we report main results in terms of community structure with the related discussion.

Dataset and main characteristics of the WTN
We refer to the World Trade Data, available on the Observatory of Economic Complexity database 3 . The database has been developed by the Research and Expertise Center on the World Economy at a high level of product disaggregation and it is based on original data provided by the United Nations Statistical Division (UN Comtrade). In particular, a harmonization procedure, that reconciles the declarations of exporters and importers, enables to extend considerably the number of countries for which trade data are available, as compared to the original dataset. In this analysis, we refer to the last version published in 2017, based on the Harmonized Commodity Description and Coding System, and that provides aggregated bilateral values of exports for each couple of origin and destination countries, expressed in billion dollars. We focus on the aggregated data of last available year, namely, 2016. Hence, we construct a directed and weighted network (see figure 5), where each node is a country and weighted links represent the amount of product trades between couple of countries.
The network is characterized by 221 nodes and 26197 links. Its density is approximatively 0.539: on average, each country has trades with more than a half of the entire network. The system is strongly connected. However, the network is not regular and is far from being complete or, in other words, most countries do not trade with all the others, but they rather select their partners. Furthermore, main trade flows tend to be concentrated in a specific sub-group of countries and a small percentage of the total number of flows accounts for a disproportionally large share of world trade. For instance, the top 10 countries export more than 50% of the total flow. The maximum weight corresponds to the channel from China to USA and its value amounts to 436 billion dollars. Minimum, non null, weights are involved in the trade between a number of very small countries, far from each others, and they are approximatively around 1 thousand dollars. A first insight in the structure of the WTN is given by the total strengths, in and out, of the nodes, which are the natural extension of the degree centrality to the case of a weighted and directed network. In this context, the in-strength s in i measures the total trade flows incoming to the country i, that is the import. The out-strength s out i measures the total trade flows outgoing from the country i, that is the export. They can be calculated by formulas 1 and 2. We observe a strict relation between in and out strength distribution with a Spearman correlation coefficient equal to 0.956. Hence, countries are ranked in a very similar way in terms of in and out strength. However, it is worth pointing out the presence of specific countries with a balance of trade significantly different from zero (see, for instance, GBR and KOR in Table 1 Table 1 Rankings of top ten countries in terms of strength (in, out and total, respectively.)

Summary of the methodology
We summarize the main steps of the methodology we are proposing. Main steps are listed in the following pseudo-code. The code has been written taking into account the communicability distance matrix Ξ, but the same procedure can be easily applied by considering the resistance matrix Ω.
1. let G be the original directed weighted network with n nodes and with a weighted adjacency matrix W; 2. symmetrise and build the undirected weighted network G 1 with weighted adjacency matrix: W 1 = 1 2 (W + W T ); 3. normalise and build the undirected weighted network G 2 with weighted adjacency matrix: W 2 = S −1/2 W 1 S −1/2 ; 4. construct the distance matrix Ξ based on the communicability distance; 5. define the threshold interval [ξ min , ξ max ], where ξ min and ξ max represent the minimum and the maximum communicability distances between couples of nodes, respectively. Set ξ h = ξ min , with h = 0. 6. define a n × n matrix M h = [m ij ] such that m ij = 1 if ξ ij ≤ ξ h and i = j 0 otherwise ; 7. build the undirected unweighted network G 3,h from the binary adjacency matrix M h ; 8. select the partition P h given by the components of the network G 3,h ; 9. define V h the set of isolated nodes in G 3,h and set G 2,h = G 2 ; 10. if the set V h is not empty, execute the following steps (a-f), otherwise skip at step 11; a) let H 2,h be the subgraph of G 2 induced by nodes that belong to V h ; b) delete from both G 2,h and G 3,h the nodes belonging to V h and all their connections; c) add the subgraph H 2,h to G 2,h ; d) set all the non-null weights of the graph H 2,h equal to 1; e) add the subgraph H 2,h to G 3,h ; f) select the partition P h given by the components of the new network G 3,h ; 11. compute the modularity of the network G 2,h with respect to the partition P h ; 12. set the number of iterations r, compute k = ξmax−ξmin r , set ξ h = ξ h−1 + k and h = h + 1 and repeat steps 6-11 until ξ h ≤ ξ max ; 13. select the optimal partition P h as the partition P h that provides the maximum modularity.
We highlight some key points of the methodology we are proposing. We aim at clustering countries on the basis of a specific distance. Varying the threshold we can disentangle the role of very tight relationships between couples of countries.
However, reducing the threshold distance a great number of isolated nodes may appear. They are typically very small countries whose trade volume is very low and whose commercial partners are few. They play a marginal role in the WTN and they do not affect in a significant way the structure of the network in terms of relevant communities. Nonetheless they contribute a lot to determine the modularity, since any isolated node in a partition tends to reduce its value. A great number of isolated nodes together with an almost completely connected giant component produce immediately negative values of modularity, preventing us from applying the maximum modularity criterion. In order to fix this problem, we propose to consider all the isolated nodes as members of a single distinct community. On one hand, this choice does not affect significantly the analysis of the core of the WTN and, on the other hand, it allows to apply the maximum modularity criterion. The treatment of isolated nodes regards step a-f of the procedure summarized before. In particular, we first identify isolated nodes in the binary network G 3,h . We extract from the original network G 2 the subgraph H 2,h induced by nodes that are isolated in G 3,h , In this way, we preserve the connections between these nodes observed in the original network. We build then two fictitious networks G 2,h and G 3,h where isolated nodes and their connections are replaced by the subgraph H 2,h . Since G 3,h is an unweighted network, we set all positive weights equal to 1. In this way we can compare the network G 2,h with a partition, extract from G 3,h , in which the isolated nodes are already well identified and collected into a single cluster (or a few clusters if the case).

Results
We initially applied the methodology described in Section 6.2 considering the communicability distance. The rationale for using the communicability metric on the WTN is the following. Two countries share a total volume of trade because they exchange a given set of products, of any kind. But they can be linked even if they don't exchange each other a given product, that is there is no direct flow of such product between them. A higher order exchange may occur between them. For instance, a country A exports some raw materials -let's say, iron -to a country B; country B produces mechanical parts from iron and exports them to country C. A and C communicate via a higher order walk and they depend on each other even if the two countries are not neighbours in the network. Communicability takes into account precisely all possible weighted walks between two nodes. Therefore, we calculate the communicability matrix G on the normalised network G 2 and the corresponding communicability distance matrix Ξ. Using this metric, we find that the nearest countries are USA and Canada with a distance ξ min = 1.242 and the farthest countries are USA and Seychelles Islands with a distance ξ max = 1.470. For each value of the threshold distance between minimum and maximum, we look at the corresponding partition in communities. In Figure 6, we plot the value of modularity (in red) and the number of communities (in blue), counting each isolated node as an independent one. Both values are expressed as functions of the threshold ξ h . The maximum of modularity (equal to 0.27) is reached at a threshold distance ξ h = 1.38. It corresponds to 21 communities, where the last one collects all the isolated nodes. We display in Figure 7 countries constituent the optimal partition and we list in Table 2 the five top biggest communities in terms of number of countries (excluding the one that groups isolated nodes).
Going deeper into the composition of the communities, community 2 (in blue) includes almost all continental European countries, with Great Britain and Ireland. This community acts on the screen of the global network as single player. It is worth pointing out also the presence of Morocco, confirming positive effects of bilateral trade agreements (see, e.g., [58]). As opposed to this community, we see the biggest community 1 (in red) which sees United States and China as main actors. This means that in Europe there are preferential channels of internal exchanges, whereas, outside Europe, most communication channels seem to be polarized around the exchange channel between China and the US and all their satellites countries. Moreover, in Europe we can recognize two other well-identified and coherent communities: Scandinavian and Baltic countries (community 3 in cyan) and the countries of the former Yugoslavia (community 4 in green). Finally, an independent role is acknowledged to the community of the countries of the Russian federation (community 5, in olive). Although a positive trade balance and a priority of Russian government of an increasing participation in the economic relations of Asia-pacific region (see [42]), at moment, results show preferential channels with border countries.  If we reduce the threshold, we let very strong channels of communication between countries emerge. For instance, Figures 8 and 9 show the community structure lowering the threshold distance (equal to ξ 0 = 1.36 and ξ 0 = 1.34, respectively). Moving from 1.38 to 1.36 some loose connections are lost. Great Britain and Ireland split up from community 2 creating a separate cluster together. The South East Asian (in violet) and former Yugoslavia (in green) disintegrate, Australia goes out from community 1, and the strong community in the South of Africa loses some country. Reducing further the threshold to 1.34, only the most closely interrelated communities survive. China disconnects from USA, but North America (in red) still behaves as a single cluster. In Europe only few inseparable relations are saved, like Spain and Portugal or Belgium, Luxembourg and Holland (in blue). Other strong links are revealed, for instance, between Brazil and Argentina, between China and Japan or between Israel and Palestine. A significant feature of our approach is the fact that it allows to get deeper insight into the internal structure of each community and to give a measure of the mutual relationships between communities. Let us refer now to the clusters depicted in Figure 7 and detected with the maximum modularity criterion. In this regard, we display in Figure 10 the distributions of the communicability distances between pair of countries that belong to the same community. In particular, we compare the distributions for the five relevant communities listed in Table 2 If we focus, for instance, on community 1 and 2, we can inspect and compare their internal structure by providing some synthetic indicators in Table  3. From the analysis of the values shown in Table 3, we can say that the community 2 (Europe) is more compact than the community 1 (USA-China) since in the former the average intracluster distance is slightly lower than in the latter. Although community 2 shows a bigger minimum intracluster distance, countries appear to be somehow more homogeneous than countries in community 1. This is partially related to the geographical distribution of the countries inside the communities. Inside the community 2, countries are all lo- cated in Europe except Morocco, which remains incorporated in the European community due to its strong commercial connection with the Mediterranean countries also thanks to the preferential agreements in force (see, e.g., [26]). Last column of Table 3 provides the same indicators computed on intercluster basis. This analysis allows to provide additional information in terms of heterogeneity in the group and between groups. Main results show high similarity between the two top communities (the average intracluster distances are 1.41 and 1.40, respectively). It is worth pointing out the lower intercluster standard deviation. It means that couple of countries that belong to a different community has a similar distance between them.  Table 3 Intercluster and Intracluster characteristics of the distributions of communicability distances. Columns Community 1 and Community 2 refer to the intracluster properties of the two main detected communities, in terms of number of components. Last column reports the corresponding intercluster properties computed between the same two communities.
It is noteworthy that additional insights can be provided by assessing the relevance of each country in the community. Indeed, communicability distance matrix provides a metric on the network and on each subnetwork, like a com- Fig. 10 Distributions of Communicability Distances between countries of the some community. We display only the distributions related to the five main communities summarized in Table 2 munity. Therefore, we adapt the idea of closeness to our context, by providing the following communicability closeness to assess how effectively a node is supposed to spread trade flows through the network. Similarly to the definition of closeness, we define the communicability closeness as: where the sum is over all the internal nodes of the cluster C to which the node i belongs.
To this end, we rank in Figure 11 countries of community 1 and community 2, respectively, on the basis of values of C i . Concerning community 1, it is worth to stress that the centre of this community is located in China, Japan and South Korea and not in the North American sub-community. The three Asian nations are nowadays major traders and their high-level economic cooperation has been strengthened also because of the speed-up of the negotiations on the trilateral Free Trade Agreement. The three parties unanimously agreed to further increase the level of trade and investment liberalization based on the consensus reached in the Regional Comprehensive Economic Partnership Agreement 4 . According to cluster 2, it is noticeable the central role that Germany plays in Europe. France has the second position and other countries like Great Britain, Italy, Spain, Benelux and Netherlands show similar relevance. Most central countries are the main economies of the European Union. As shown in [13], these countries are the core group of the network. Moreover, it is interesting to see that most central country in a community has not necessarily the same relevance on the whole network. We have indeed that, in terms of subgraph centrality, when we deal with the whole network (see Figure 12), USA appears as the key player followed by China and Germany. This ranking is inline with the top three provided by the World Trade Organizations [64] in terms of World's leading traders of goods and services [64].
Additionally, it is interesting to highlight that the relevance of countries reported in Figure 12 is consistent with the Economic Complexity Index (ECI), introduced by [38]. The ECI allows to rank countries in the WTN according to the diversification of their export flows, which reflects the amount of knowledge that drives their growth. The higher is the ECI, the more advanced and diversified is an economy. In particular, countries whose economic complexity is greater than expected (on the basis of their global income), tend to grow faster than rich countries with a low ECI. In this perspective, ECI represents a suitable tool for comparing countries in the WTN independently of their total output and it has been extensively validated as a relevant economic measure by showing its capability to predict future economic changes and to explain international differences in countries incomes. Although the network we analysed in the present work is based on the total output and this fact prevents us from comparing directly their values with the ECI for a given country, there is a positive correlation between them. All the top 20 countries in Figure 12 show a positive and high value of ECI. More specifically, they kept a high value of ECI during the years preceding the year to which the network refers (2016) and this can justify the high value in the aforementioned centrality measures.
Finally, from the point of view a single country, it is worth to look for the closest trade partners, that is the nearest nodes in terms of communicability distance. Figures 13 show the distance profiles for the most central country of community 1 and community 2, respectively. For instance, looking at Figure 13 (right-hand side), we can notice countries, as Austria, Poland, Czech Republic that are characterized by a condition of strong dependence on Germany, that is a major player in the network. Similarly Figure 13, left-hand side, shows how strong is the commercial relationship between China and Hong Kong, also as a result of the trade agreements between the two countries, like CEPA (Closer Economic Partnership Arrangement) aimed at eliminating duties on large categories of products. For the Chinese trade market Hong Kong plays a crucial role since foreign companies use Hong Kong as a springboard to invest in China thanks to its infrastructure network that has no equal in the world, investor protection, transparent and efficient judicial system, legal certainty.

Results in terms of resistance metric
The methodology described in Section 6.2 has also been applied using the resistance distance ω. In this case, we consider the total trade of a given country as flows of the global wealth that has been produced during a year. Therefore, the Gross Domestic Product (GDP) is the attribute of interest on each node. In this regard, the effective resistance of an edge expresses how easily (or not) a unit flow moves from a country to another one, i.e. how easily two countries trade a unit of wealth, independently of its nature. It is noteworthy that, according to formula 10, the resistance distance between a pair of countries depends on the values of the vibrational centralities of both countries and on the value of their mutual connections.
Therefore, we construct the vibrational communicability matrix G v on the normalised network G 2 , and the corresponding resistance distance matrix Ω. Using this metric, we find that the nearest countries are, again, USA and Canada with a distance ω min = 1.53 and the farthest countries are USA and Germany with a distance ω max = 2.24. For each value of the threshold distance between minimum and maximum, we obtain the corresponding partition in communities. The maximum modularity partition corresponds to 16 final communities, where, as before, isolated nodes have been grouped all together in the same community. The value of modularity for this partition is Q = 0.31. In Figure 14, we plot the value of modularity in red and the number of communities, counting each isolated node as an independent one, in blue as functions of the threshold ω 0 . The maximum of modularity is reached at a threshold distance ω = 1.85. Main results in terms of geographical distribution are displayed in Figure 15 and, as in the previous Section, we summarize in Table 4 main composition of top communities in terms of number of constituents.
With respect to results based on communicability, we have that the first community has a larger number of countries (equal to 66). Additionally, the larger community includes again main asian and oceanian countries but North America behaves as a separate cluster. This result is in line with the literature that emphasizes the interesting economic relation between Asia and Oceania. Several works showed that the Asia-Oceania community collapsed after China entered the WTO in 2001 and built strong trade relationships with other communities, especially with the external cores, (i.e. the United States and Germany). China then became regionally attractive and restored the Asia-Oceania community as the community leader after it gained a significant portion of trade globally (see, e.g., [65]). Similar results to the communicability case are observed for the main European community. However in this case, Great Britain and Ireland provide a separate group and Italy and Morocco are isolated countries.
It is worth pointing out that communities detected above represent groups of countries showing a positive correlation in their trade strength, whereas members of different clusters show a negative correlation. Being strongly anticorrelated means that when the total trade deficit of a country grows, the total trade surplus of a second country grows too. For instance, Japan and USA have been classified in different communities by the methodology. Indeed, in the literature, empirical analyses show a negative correlation coefficient between normalised trade strengths of these countries (see, e.g., [41] and [51]). Similar arguments can be extended also to other pairs of countries. For instance, Germany is negatively correlated with USA (see [41]) and show a high positive correlation with Belgium and France (see [51]), that belong to the same community.
If we disentangle communities characterized by very tight relationship between countries, the results seem strictly related to the ECI index. We may expect that, if two countries communicate well, then their ECI's could be similar. That is, if their mutual distance is small, both in terms of communicability metric and resistance metric, then they display similar values of ECI. In fact, the existence of multiple channels of trade exchange between them would result in a similar diversification of their output. This means that countries inside each community (could) share homogeneous values of ECI. Concerning Table  4, we notice small clusters whose components show homogeneous values of the ECI index.  As in the previous Section, we explore main characteristics of two most relevant communities (see Table 5), It is noticeable that, although the two groups show a very similar distance, European countries are characterized by a higher heterogeneity. Focusing on intercluster indicators, we notice a lower  Table 4 Members for the six main communities in terms of number of countries. The optimal partition has been obtained by applying the procedure based on the maximum modularity and a threshold depending on the resistance distance.
similarity between the two communities with respect to Table 3 based on communicability.  Table 5 Intercluster and Intracluster characteristics of the distributions of resistance distances. Columns Community 1 and Community 2 refer to the intracluster properties of the two main detected communities, in terms of number of components. Last column reports the corresponding intercluster properties computed between the same two communities.
The relevance of the country can be now assessed in terms of vibrational centrality. To this end, we display in Figure 16, the top 20 countries, calculated over the whole network. China, USA and Germany are again in the top 3, with China playing as the best spreader node. Also in this case, almost all the top 20 has a positive ECI. A comparison between Figures 12 and 16 confirms the different role played by USA and China in the global network. As confirmed by [64], USA is the leading commercial service provider and in such a way it is widespread well-integrated in the global market; on the other side, China plays the role of hub for goods and represents the leading merchandise trader and this gives to the country a very robust position which makes it less vulnerable to market turmoil.
Finally, from the point of view a single country, it is worth to look for the closest trade partners, that is the nearest nodes in terms of resistance distance. Figure 17 shows the distance profiles for the most central country of community 1 and 2, respectively. These plots can be interpreted as the list, in decreasing order, of countries that are most positively correlated with the selected centre, China or Germany. For instance, while in terms of communicability distance China is well-communicating with USA (third position in Figure 13), USA does not belong to the top 20 most correlated countries with China. Rather, the left-hand side in Figure 16 clearly shows a driving and synchronizing effect of the Chinese giant in the entire South-East Asia area. Similarly, Figure 16, right-hand side, confirms the role of Germany in the European Union and the strong correlation with Austria, Czech Republic and Poland.

Conclusions and further research
Community detection is a key topic in the analysis of complex systems, where discovering the inner structure plays a relevant role. In particular, the centrality of countries and the relationships between them assume specific relevance in the World Trade Network, where economical and geopolitical phenomena affect over time the structure of the global network. In this framework, this work aimed at detecting different levels of clustered communities in the network on the basis of both communicability and resistance distances. The proposed methodology allows to discover the hidden hierarchical structure of the network, as it presents a degree of flexibility highlighting very tight relationships by varying the threshold parameter, and revealing in this way the clusters of nodes that more easily communicate. Moreover, it performs well also for Features and properties of each community can be exploited in order to compare the characteristics of different clusters and to detect the most central countries inside the single community as well in the whole network. Numerical results depict the structure of the economic trade detecting main relevant communities. In particular, main community sees United States and China as main actors. Most flows are polarized around the exchange channel between China and USA and all their satellite countries. However, focusing on the correlation between trades, the procedure emphasizes the different role of these two countries. In particular, it is worth mentioning the emerging of China-Oceania community when deep links emerge. Furthermore, it is confirmed that Germany plays a key role in Europe and preferential channels of internal exchanges are observed in the European market. be the flow of such an attribute from node i to node j. We denote by I i the total outgoing flow from the node i to its adjacent nodes, i.e. I i = n j=1 a ij (v i − v j ). In matrix form, the total outgoing flow of the nodes attribute is then The Laplacian matrix transforms nodes attributes v i , i = 1, ..., n into outgoing flows from nodes I i , under the assumption that a flow I ij along a given edge is equal to the gradient ∆v ij = v i − v j . This assumption is equivalent to choose an effective resistance equal to 1 along all edges. Of course, we may have both outgoing and ingoing currents according to the sign of ∆v ij : positive for outgoing flows from i and negative for ingoing flows into i. A similar meaning can be given to L + 1 n J v. Indeed, where v = 1 n n k=1 v k , that is, the operator L + 1 n J adds to the flows a constant term given by the mean value of all the attributes of the nodes. Then, the matrix L + 1 n J transforms nodes attributes v into total outgoing flows I in the network, up to an additive constant. In a similar way, the inverse L + 1 n J −1 acts on a current vector I and produces a state vector v, which can be interpreted as the cause of such currents in the network. Specifically where, once again, the term 1 n JI = Iu is the average value of the outgoing currents coming from every node.
Suppose now that in the system there are an outgoing flow equal to 1 from a node (node 1, for instance), an ingoing flow equal to −1 into another node (for instance, node 2), whereas for all the other nodes the flow is zero. This is equivalent to a current vector equal to I = [1, −1, 0, . . . , 0] T = e 1 − e 2 . Loosely speaking, a unit information is coming out from node 1 and goes entirely into node 2. To produce these flows, we have to start from an initial attributes vector on nodes given by where the last equality holds because J(e 1 − e 2 ) = 0, that is I = 0. Thus, the resistance distance between nodes 1 and 2 is given by ω 12 = (e 1 − e 2 ) T L + 1 n J −1 (e 1 − e 2 ) = v 1 − v 2 = ∆v 12 .
If ∆v 12 is small, a small gradient is enough to transmit such a unit flow from node 1 to node 2; whereas, if v 1 − v 2 is big, a high gradient is needed in order to produce the same unit flow. More in general, let's imagine that in the node i the value v i is positive. Then the fact that another attribute v j with j = i is positive means that node i and node j are strongly correlated since it is enough a low attribute difference to subtract from node i a unit flow. This means that these two nodes communicate a lot. Whereas, if for another node k with k = i, the corresponding component v k is negative this implies that node i and node k are strongly anti-correlated since, in order to produce a unit flow from node i, node k has to be at a negative attribute, i.e. the attribute difference between i and k must be high. This means that the two nodes don't communicate well. The signs of the components of the vector v indicate nodes that are positively or negatively correlated with node i according to the fact these components have the same sign as v i or not. Let us observe that, in general, v = L + I = L + (e i − e j ) = L + i − L + j with L + i i-th column of the matrix L + . That is, if we want to decrease by 1 the attribute of node i and increase by 1 the attribute of node j, we have to take an initial distribution of attributes on nodes equal to the difference between i-th column of L + and j-th column of L + , and these columns are also the values of vibrational communicability G v between nodes, as defined in the text.


Solving the Asymmetry Multi-Objective Optimization Problem in PPPs under LPVR Mechanism by Bi-Level Programing

: Optimizing the cost and beneﬁt allocation among multiple players in a public-private partnership (PPP) project is recognized to be a multi-objective optimization problem (MOP). When the least present value of revenue (LPVR) mechanism is adopted in the competitive procurement of PPPs, the MOP presents asymmetry in objective levels, control variables and action orders. This paper characterizes this asymmetrical MOP in Stackelberg theory and builds a bi-level programing model to solve it in order to support the decision-making activities of both the public and private sectors in negotiation. An intuitive algorithm based on the non-dominated sorting genetic algorithm III (NSGA III) framework is designed to generate Pareto solutions that allow decision-makers to choose optimal strategies from their own criteria. The e ﬀ ectiveness of the model and algorithm is validated via a real case of a highway PPP project. The results reveal that the PPP project will be ﬁnancially infeasible without the transfer of certain amounts of exterior beneﬁts into supplementary income for the private sector. Besides, the strategy of transferring minimum exterior beneﬁts is more beneﬁcial to the public sector than to users.

Introduction
Public-private partnerships (PPP) have been applied worldwide as a financing instrument to provide public infrastructures and services. In a PPP project, the private sector is usually responsible for the operation of the project and for providing part of the initial investments, with a return on governmental subsidies, user payments or extraneous earnings derived from exterior benefits generated by the project [1]. In other words, when the private sector benefits increase, the costs from the public side are increased.
In the past few decades, the view that cost or benefit allocation among multiple players in a PPP project is actually a multi-objective optimization problem (MOP) has been formed. Koppenjan (2005) investigated the reasons for the stagnation of transport PPP in the Netherlands, finding that the mutually exclusive goals of the public and private sectors affected the sustainability of cooperation [2]. Papajohn (2011) put forward the view that the public is an important stakeholder of PPPs whose interest would be damaged if the government uses PPP to make excessive investments [3]. Garvin (2010) believed that PPP decision-making from the government needs to justify its economic feasibility by

The LPVR Mechanism
The original idea of LPVR is a flexible-term concession; in other words, the concession term is tied to certain return on investment (ROI) indicators for the private sector [17]. When future demand is better than expected, the private sector is allowed to end the contract earlier than estimated. In contrast, when the future demand is not as optimistic as expected, the private sector is allowed to extend the concession term in order to achieve the pre-set ROI. According to Engle et al. (1997,2001), the pre-set ROI will be used at a discounted rate to calculate the private sector's present value of revenue (PVR), which has been taken as the key variable in the tender of the concession contract [15,18]. The private sector with the least PVR, i.e., the LPVR, will be granted the concession contract. Equation (1) presents the general expression of the LPVR [18]: where pri LPVR is the net present value calculated in the base year; 0 C is the initial investment that is paid by the private sector in a PPP project; τ is the estimated concession term; ρ is the unit charge from users; i q is the estimated quantity of service or products provided by the project in the i-th year; and i m is the maintaining and operating cost in i-th year. Theoretically, an LPVR of 0 means that the private sector obtains revenues that are exactly equal to the costs incurred.
In the context of a PPP project, the net present value of revenue of the private sector-the PVR of the private sector-is also affected by decisions from the public sector, which cause variations in both the initial investment section and operation section. Such variations are shown in Equation (2), where (1) Δ is the initial investment shared by the public sector and (2) i Δ contains the operation subsidy and consumed exterior benefits of the i-th year. During the negotiation, the private sector needs to find the least PVR (the LPVR) to win in the bidding process.

The LPVR Mechanism
The original idea of LPVR is a flexible-term concession; in other words, the concession term is tied to certain return on investment (ROI) indicators for the private sector [17]. When future demand is better than expected, the private sector is allowed to end the contract earlier than estimated. In contrast, when the future demand is not as optimistic as expected, the private sector is allowed to extend the concession term in order to achieve the pre-set ROI. According to Engle et al. (1997,2001), the pre-set ROI will be used at a discounted rate to calculate the private sector's present value of revenue (PVR), which has been taken as the key variable in the tender of the concession contract [15,18]. The private sector with the least PVR, i.e., the LPVR, will be granted the concession contract. Equation (1) presents the general expression of the LPVR [18]: where LPVR pri is the net present value calculated in the base year; C 0 is the initial investment that is paid by the private sector in a PPP project; τ is the estimated concession term; ρ is the unit charge from users; q i is the estimated quantity of service or products provided by the project in the i-th year; and m i is the maintaining and operating cost in i-th year. Theoretically, an LPVR of 0 means that the private sector obtains revenues that are exactly equal to the costs incurred.
In the context of a PPP project, the net present value of revenue of the private sector-the PVR of the private sector-is also affected by decisions from the public sector, which cause variations in both the initial investment section and operation section. Such variations are shown in Equation (2), where ∆ (1) is the initial investment shared by the public sector and ∆ (2) i contains the operation subsidy and consumed exterior benefits of the i-th year. During the negotiation, the private sector needs to find the least PVR (the LPVR) to win in the bidding process.

The Interactive Stackelberg Decision-Making Process
Equation (2) shows that the private sector benefits are affected by the public sector's decisions. Correspondingly, the public sector benefits are also affected by the decisions of the private sector when the control variables ρ and τ of the private sector are used as input parameters for calculating benefits (or costs) at the public side. Supposing ∆ (1) and ∆ (2) i in Equation (2) are determined by strategy u proposed by the public sector, and ρ and τ are the elements of strategy v proposed by the private sector, an interactive decision-making process with the public sector as the leader and the private sector as the follower is shown in Figure 2.

The Interactive Stackelberg Decision-Making Process
Equation (2) shows that the private sector benefits are affected by the public sector's decisions. Correspondingly, the public sector benefits are also affected by the decisions of the private sector when the control variables ρ and τ of the private sector are used as input parameters for calculating benefits (or costs) at the public side. Supposing (1) Δ and (2) i Δ in Equation (2) are determined by strategy u proposed by the public sector, and ρ and τ are the elements of strategy v proposed by the private sector, an interactive decision-making process with the public sector as the leader and the private sector as the follower is shown in Figure 2. The above interactive decision-making process could be characterized by the Stackelberg decision-making problem [19], which has the following characteristics in general: (1) each decisionmaker has specific objectives and control variables, and their strategy will affect the benefits of the others as well as their own; (2) the objectives of different decision-makers are located in different levels, presenting a hierarchical decision-making structure-the upper decision-maker has higher priority in terms of the action order than the lower decision-maker; (3) the decisions of the upper and lower level impact each other, thus forming a relationship of mutual influence and restriction; and (4) the final decision will be a globally feasible solution that satisfies all decision-makers as much as possible.
The Stackelberg decision-making problem could be solved by a bi-level programing (BLP) model [20]. The general form of the BLP model consists of an upper programing (UP) problem and lower programing (LP) problem, which could be expressed as Equations (3)-(6), as shown in [21].
(UP) (5) subject to The above interactive decision-making process could be characterized by the Stackelberg decision-making problem [19], which has the following characteristics in general: (1) each decision-maker has specific objectives and control variables, and their strategy will affect the benefits of the others as well as their own; (2) the objectives of different decision-makers are located in different levels, presenting a hierarchical decision-making structure-the upper decision-maker has higher priority in terms of the action order than the lower decision-maker; (3) the decisions of the upper and lower level impact each other, thus forming a relationship of mutual influence and restriction; and (4) the final decision will be a globally feasible solution that satisfies all decision-makers as much as possible.
The Stackelberg decision-making problem could be solved by a bi-level programing (BLP) model [20]. The general form of the BLP model consists of an upper programing (UP) problem and lower programing (LP) problem, which could be expressed as Equations (3)-(6), as shown in [21].
(UP) min subject to g(u, v) ≤ 0 (6) where F(u, v(u)), u and G(u, v(u)) are the objective function, vector of control variables and constraint of the upper level decision-maker, respectively; f (u, v), v and g(u, v) are those of the lower decision-maker; and u is used as the input parameter when calculating f (u, v) and v(u) is used as the input parameter when calculating F(u, v(u)).

Modeling the Problem by Bi-Level Programing (BLP)
There are two prerequisites for both parties to participate in the PPP project: one is that the public sector should accomplish "value for money" (VFM)-the lifecycle costs paid in the PPP mode should be no more than that in the traditional investment mode (without private funding)-while the other one is that the private sector can obtain an actual return level that is no less than the pre-set rate of return. In other words, the values of VFM and LPVR are required to be positive to satisfy the two prerequisites.
Definitions of parameters and variables in the formulas of the model are listed in Table A1 of Appendix A.

The Upper Programing (UP) Problem
The UP problem is a three-objective MOP that is used to minimize the lifecycle cost of the public sector, the payment level of users and the consumption of exterior benefits. The constraint for the public sector to participate in a PPP project is to realize "value for money", i.e., the value of VFM should be non-negative. The objectives and constraints of the UP problem are formally expressed as follows: where G α, δ, ν (α,δ,ε) (ρ, τ) is the net present value of the lifecycle cost paid by the public sector under the PPP mode; ρ is the unit toll per vehicle-kilometer, representing the payment level of users; and ε is the ratio of collected exterior benefits used as supplementary income for the project. ν (α,δ,ε) (ρ, τ) is the input variable combination of ρ * , τ * , which is solved by the lower programing (LP) model upon the upper decision {α, δ, ε}. The full expression of the UP objective is: which consists of two terms. The first term of Equation (11) is the net present value (NPV) of equity funds shared by the public sector, and the second term is the NPV of subsidies that the public sector expends during the whole operation period. The value of VFM is defined as the difference between the lifecycle cost under traditional investment mode, called the public sector comparator (PSC) value, and that under the PPP mode expressed by G α, δ, ν (α,δ,ε) (ρ, τ) . Under the traditional investment mode, the equity funds of construction investment and operation costs are all afforded by the public sector. The value of PSC can be expressed as: Ir cap r con (1−α)Ir cap r con

The Lower Programing (LP) Problem
In the LP problem, the private sector company reduces their benefits due not to an original intention but to improve the chance of being selected as the partner, especially in a competitive bidding scenario under the LPVR mechanism. The objective function of the private sector company is to minimize their PVR, which is constrained by their own non-negative requirement. The model of the LP problem is formally expressed as follows: where P(α , δ , ε , ρ, τ) is the net present value of the private sector benefit, and {α , δ , ε } is given by the upper decision-maker (the public sector). The full expression of P(α , δ , ε , ρ, τ) is:

Optimization of Objectives in the BLP Model
As is known, it is difficult to find an absolute optimal solution in a MOP with conflicted objectives in which a series of comparatively optimal solutions, called Pareto solutions, could be obtained. Each Pareto solution performs better for at least one objective and worse for at least one objective than other Pareto solutions [22]. Pareto solutions can leave the decision-makers with adequate space for strategies without any requirements of weighting the multiple objectives as one general objective in advance.
In this paper, Pareto solutions will be solved and then compared to the original case in a case study to determine which objectives are present and the degree to which they are optimized. Figure 3 shows the schematic of how the Pareto solutions are demonstrated in the final results. For instance, in Figure 3, Pareto solution 1 reduces the payment of the public sector and users but increases the consumption of exterior benefits (EB) and the LPVR of the private sector; Pareto solution 2 reduces the payment of the public sector and users, keeps the consumption of EB the same, and lifts the LPVR; and Pareto solution 3 reduces the public sector costs, keeps the user payment level the same, and increases the consumption of EB and LPVR. It is noticeable that when the value of LPVR is more than zero, this means that the private sector could obtain more revenue than expected.

NSGAIII Framework-Based Algorithm for the BLP Model
BLP models are difficult to solve directly because of their status as NP-hard problems. The model built in Section 3 adds even difficulty since the UP problem is a three-objective nonlinear MOP and the LP problem is a mixed-integer nonlinear problem. Interactions between the two levels are hardly predicted; thus, the widely used Karush-Kuhn-Tucker (KKT) conditions [23] or sensitivity analysisbased methods [24] would not work for this problem.
In studies [25][26][27][28] aimed at solving such complex BLP models in recent years, embeddedcomputational approaches are applied. In general, the outer layer of their framework contains a heuristic algorithm to solve the UP problem, such as the genetic algorithm (GA), simulated annealing (SA) and particle swarm optimization (PSO). Then, a specific algorithm is applied to provide the optimal solutions to the LP problem, which are used as input parameters to compute the UP problem.
In this paper, the general idea used to solve the complete BLP model is to create an embeddedcomputational framework. The algorithm of non-dominating GAs-known as NSGAIII-is used to solve the three-dimensional MOP in the upper level. Additionally, the algorithm based on variable discretization is taken to quickly produce optimal solutions to the LP problem. NSGAIII is an algorithm developed on the GA framework to solve many-objective problems (MaOPs) [29] and has a good reputation for improving the evenness and diversity of solutions via the elite selection strategy of "reference points" [30].
In the built BLP model, problems α δ ε is an optimal solution of the upper problem. There are two main targets of the algorithm design of the BLP model: (1) to ensure that the LP problem produces feasible solutions upon each decision from the UP; and (2) to create a fast searching technique, which is required to solve the LP problem to avoid embedded iterations that would consume too much computing time, since the algorithm of the UP problem is evolution-based.

Interception Strategy for Reparing Solutions of the UP Problem
The main idea of this strategy is to intercept the UP solution{ } , , α δ ε that causes the LP not to be able to find a feasible solution under its own constraints. The interception strategy process includes the following steps.

NSGAIII Framework-Based Algorithm for the BLP Model
BLP models are difficult to solve directly because of their status as NP-hard problems. The model built in Section 3 adds even difficulty since the UP problem is a three-objective nonlinear MOP and the LP problem is a mixed-integer nonlinear problem. Interactions between the two levels are hardly predicted; thus, the widely used Karush-Kuhn-Tucker (KKT) conditions [23] or sensitivity analysis-based methods [24] would not work for this problem.
In studies [25][26][27][28] aimed at solving such complex BLP models in recent years, embedded-computational approaches are applied. In general, the outer layer of their framework contains a heuristic algorithm to solve the UP problem, such as the genetic algorithm (GA), simulated annealing (SA) and particle swarm optimization (PSO). Then, a specific algorithm is applied to provide the optimal solutions to the LP problem, which are used as input parameters to compute the UP problem.
In this paper, the general idea used to solve the complete BLP model is to create an embedded-computational framework. The algorithm of non-dominating GAs-known as NSGAIII-is used to solve the three-dimensional MOP in the upper level. Additionally, the algorithm based on variable discretization is taken to quickly produce optimal solutions to the LP problem. NSGAIII is an algorithm developed on the GA framework to solve many-objective problems (MaOPs) [29] and has a good reputation for improving the evenness and diversity of solutions via the elite selection strategy of "reference points" [30].
In the built BLP model, problems of both levels have independent constraints. Every individual in the final Pareto solution set of α, δ, ε, ρ, τ should be globally feasible, i.e., satisfying the constraints of both the UP problem and the LP problem. In each part of the combination of α, δ, ε, ρ, τ , ρ, τ is solved based on {α, δ, ε}, and {α, δ, ε} is an optimal solution of the upper problem. There are two main targets of the algorithm design of the BLP model: (1) to ensure that the LP problem produces feasible solutions upon each decision from the UP; and (2) to create a fast searching technique, which is required to solve the LP problem to avoid embedded iterations that would consume too much computing time, since the algorithm of the UP problem is evolution-based.

Interception Strategy for Reparing Solutions of the UP Problem
The main idea of this strategy is to intercept the UP solution {α, δ, ε} that causes the LP not to be able to find a feasible solution under its own constraints. The interception strategy process includes the following steps. Step 1: Generate N individuals of α, δ, ε, ρ max , τ max by the GA operator. ρ max and τ max are the upper limits of ρ and τ, respectively. N is the size of the population.
Step 2: Compute the value of the LP constraint function with N individuals of α, δ, ε, ρ max , τ max . Count the number of individuals which violate the LP constraint, and notate it as n.

Subsection Approximation Strategy for Optimizing the LP
When solving the LP, {α, δ, ε} is taken as the given parameter set. Besides this, τ is a set with a limited number of integer elements indicating the year of concession term, as τ= {1, 2, · · · , τ max }. Thus, Equation (16) can be equivalently transformed into Equation (17): Figure 4 shows the three-dimensional surface diagram with contours of the functional relationships between F LP , ρ and τ, where ρ is taken as a discretized variable. The F LP value remains the same at different positions on the same contour line. It can be seen from Figure 2 that larger values of ρ and τ produce larger values of F LP . Besides this, under the same value of F LP , ρ and τ present a kind of inverse correlation. Assuming that the private sector is more averse to long-term concessions than charging lower prices, the process of the subsection approximation strategy includes the following steps.
Symmetry 2020, 12, x FOR PEER REVIEW 8 of 19 Step  Thus, Equation (16) can be equivalently transformed into Equation (17): Figure 4 shows the three-dimensional surface diagram with contours of the functional relationships between LP F , ρ and τ , where ρ is taken as a discretized variable. The LP , ρ and τ present a kind of inverse correlation. Assuming that the private sector is more averse to long-term concessions than charging lower prices, the process of the subsection approximation strategy includes the following steps. Step 1: Discretize the continuous decision variable Step 2: According to the order that  Step 1: Discretize the continuous decision variable ρ ∈ [ρ min , ρ max ] according to some acceptable precision. For example, ρ min = 0 and ρ max = 0.69 at a precision of 0.01; then, ρ can be segmented into 70 pieces as ρ = {0.00, 0.01, 0.02, · · · , 0.69}.
Step 4: The finally obtained α, δ, ε, p , τ is taken as the approximate optimal solution of the LP.

Algorithm Design for the Complete BLP Model
The entire process of the proposed algorithm for the built BLP model is implemented with the 'PlatEMO' [31] source code, including the following steps ( Figure 5):   Step 0 Initialize the parameters of the NSGAIII algorithm. These parameters are population size N, the number of objective functions M and the maximum generation times gen. The actual size of population N will be determined by the inherent reference point generation method of NSGAIII [32].
Step 1: Generate a population of N individuals of {α, δ, ε} by real coding. In other words, generate N random numbers between the upper limit and lower limit of each variable (α, δ or ε), then rearrange them as an N × 3 matrix. Then, perform the interception strategy listed in Section 4.1 to obtain N globally feasible solutions, where each of them is a combination of α, δ, ε, p, τ .
Step 2: Calculate values of the UP problem's objective functions with N α, δ, ε, p, τ and find the minimum value for each objective to form the coordinate set Z min of the "ideal point". For example, each of z Step 3: Generate N parent individuals of {α, δ, ε} by the tournament selection strategy [33], then perform the Simulate Binary Crossover (SBX) [34] and Polynomial Mutation (PM) [35] operations to obtain the child population for the current epoch. Repeat the strategy listed in Section 4.1 with the tournament, SBX and PM operations as needed to make all child individuals available to solve the solutions of the LP.
Step 4: Combine the parent population and child population together to obtain 2N individuals of α, δ, ε, p, τ ; then, calculate their objective functions as fitness.
Step 5: Perform the Effective Non-Dominating Sorting (ENS) strategy [36] on the 2N individuals and allocate a non-dominating serial number to each individual. The elite selection strategy based on reference points [37] should be carried out for those individuals with the largest serial number.
Step 6: Repeat Step 3 to Step 5 until the maximum generation time gen is reached. The final population of individuals is taken as the Pareto solution set of the entire BLP.

The Big Outer Ring Highway PPP Project
According to the PPP implementation plan of the Big Outer Ring Highway (Figure 6), the total investment (I) is about 13.5 billion Yuan (1.94 billion USD at an exchange rate of 1:6.9501), the construction period is 2 years (2017-2018) and the operation period is 25 years (2019-2043). The capital fund ratio (r cap ) is 25% of the total investment, and the project has been granted a sum of extra funding (I ) from the central government equal to the amount of the capital funding. The interest rate of long-term loans for more than five years (r deb ) is 4.9%. The expected ratio of non-toll income to toll income (r nt ) is 2%. The least expected ROI of the private sector (r pri ) is 8%, and the ROI of the public sector (r gov ) is 4.25%. Other input parameters of the operation period are listed in Table 1.
The upper limit and lower limit of the decision variables are shown in Table 2. The upper limit of variable δ is determined according to the maximum unit traffic subsidy (shown in Table 2, column 4), while other values are extracted from the project plan report.
1 Figure 6. Sketch of the location of the Big Outer Ring Highways on a plain map.   Table 2. Values of the upper limit and lower limit of decision variables.

Concession Term τ (Years)
Upper limit

Evaluation Effects of Different Population Sizes
Population size is one of the key parameters of the non-dominating sorting GA algorithms. The effects of two population sizes N = 60 and N = 120 are compared with the NSGAIII algorithm. Figure 7 shows the distribution of the Pareto solutions of the two population sizes from different views in the three-dimensional space. Arrows point to the convergence direction. It can be seen from Figure 7 that the group of N = 120 performs better in terms of the convergence and evenness of solutions. There are also some solutions of group N = 60 that are relatively far away from the front edge (in the convergence direction) of the solution set, which means that these solutions are inferior to those that are right on or close to the front edge.
It can be seen from Figure 8a that a larger population creates better convergence, even though there is little difference between the two sizes of population in general. Values of the LP objective functions in the first two iterations are much higher than those in the subsequent iterations; thus, the convergence curve starts from the third to the 50th iteration, as shown independently in Figure 8b. It can be seen in Figure 8b that a larger population size leads to a better convergence than a smaller population size after the 17th iteration. Results of the group are kept for analysis in later sections of this paper.

Results and Discussion
gov C , user C , ε , VFM, LPVR and the private sector's IRR for all the 120 Pareto solutions, arranged in the order of increasing LPVR values, are plotted as stacked lines in Figure 9a. From the 120 Pareto solutions, the optimal strategies for minimizing the costs of the public sector (S1), minimizing user payment (S2), minimizing the transferred exterior benefits (S3) and minimizing the LPVR of the private sector (S4) are drawn out separately in Figure 9b. For each type of optimal strategy, only the strategies which generated the lowest LPVR are chosen. In that case, all four optimal strategies are located within the range of the first 10 solutions.

Results and Discussion
C gov , C user , ε, VFM, LPVR and the private sector's IRR for all the 120 Pareto solutions, arranged in the order of increasing LPVR values, are plotted as stacked lines in Figure 9a. From the 120 Pareto solutions, the optimal strategies for minimizing the costs of the public sector (S1), minimizing user payment (S2), minimizing the transferred exterior benefits (S3) and minimizing the LPVR of the private sector (S4) are drawn out separately in Figure 9b. For each type of optimal strategy, only the strategies which generated the lowest LPVR are chosen. In that case, all four optimal strategies are located within the range of the first 10 solutions.  It can be seen from Figure 9a that VFM and LPVR are both less than zero in the original case (the red dot lines), illustrating that strategies of the original case are not feasible for implementing the project via PPP. The ratio of consumed exterior benefits in all Pareto solutions is above that in the original case, which means the problem might be infeasible without transferring certain amounts of exterior benefits as supplementary income into the project. Variation of C gov and of VFM are not linearly correlated, revealing their sum, i.e., the PSC value is not stable. Therefore, it is not appropriate to use the VFM value to evaluate different PPP strategies.
Tables 3-5 separately show the changes when each optimal solution is compared with the original case solution (S 0 ) in different indicators. Their definitions are listed in Table A2 of Appendix A. ∆ 1~∆4 show the increments of S 1~S4 when compared to S 0 . Changes in α, ε, irr pri and τ are measured in absolute increments as: Changes in δ, ρ, C gov , PSC, VFM, F equity gov , S r gov , and B pri ,R sub ,R tol , F equity pri and C ope are measured in relative increments as: The following discussion is based on the results shown in Tables 3-5. S 1 offers a scheme of nearly complete privatization, which is without any capital fund or subsidy paid by the public sector. The results of S 1 reveal that the public sector can achieve maximum benefit when delivering the project in a fully market-based schema.
S 2 provides a schema equivalent to opening up the project to the public for free. Besides, adding 29.43% of exterior benefits realizes a 46.5% reduction in the public sector cost.
S 3 shows that at least 3.55% of exterior benefits need to be transferred into the project to meet the two prerequisites of PPP. In S3, the unit payment of users increases by 1.75% while the public sector cost reduces by 14.75%, which indicates that the strategy of minimum consumption of exterior benefits is more beneficial to the public sector than to users.
In S 4 , the actual IRR of the private sector is almost the same as the other three optimal strategies. Compared with S 1 , S 2 and S 3 , S 4 provides a more balanced schema, in which the consumption of exterior benefits increases by 30.93%, the cost of the public sector reduces by 66.62% and the unit payment of users decreases by 69.03%.
It can be seen in all four optimal strategies that shortening the concession term plays the most critical role in cutting down the cost of the public sector and increasing the benefit of the private sector. Both parties benefit from the reduction of their own costs in the operation period.

Conclusions
This paper studies the problem of Stackelberg game optimization under an LPVR mechanism when negotiating a transport PPP project between the public sector and the private sector. The decision-making scenario and the characteristics of the Stackelberg game during procurement are described and analyzed. Based on that analysis, a BLP model is constructed. The objective of the UP problem is to minimize the costs, including the public sector cost and the user payment. The constraint of the UP is to realize value for money for the public sector. The objective of the LP problem is to minimize the difference between the actual returns and the expected returns of the private sector, and its constraint is that the actual returns must be no less than the expected level.
An intuitive algorithm based on the NSGAIII framework is designed to solve the BLP model. In order to generate global feasible solutions that satisfy both the UP and LP constraints, an interception strategy is carried out to avoid infeasible solutions from UP to LP. Based on the analysis of the objective function of the LP, a subsection approximating algorithm to solve the LP by discretizing variables of the concession term and toll is proposed.
The effectiveness of the model and algorithm is validated by the PPP project of the Beijing Big Outer Ring highway. Solutions maximizing each single objective are selected from the Pareto solution set. Costs and benefits for both the public sector and the private sector under different solutions are compared and discussed. Main findings include that (1) the PPP project will be financially infeasible without transferring certain amounts of exterior benefits into supplementary income for the private sector; (2) the strategy of transferring the minimum exterior benefits is more beneficial to the public sector than to users; (3) the public sector can achieve maximum benefit when delivering the project in a fully market-based schema; and (4) the value of VFM is not suitable for use as the evaluation criteria for PPP investment schemas since its computation base, the PSC, is not static.
It is necessary to know that results in this study are case-specific, as it has not been ensured that each project generates the same level of returns. This paper is models the Stackelberg decision-making problem of two parties only. In future research, the public, local residents, tax payers, pure investors and other stakeholders in PPPs could be involved, and a multi-layered Stackelberg MOP could be built to further simulate more complex decision-making scenarios.
Author Contributions: F.L. drafted the original manuscript; J.L. and X.Y. supervised the entire study work and writing of the paper. All authors have read and agreed to the published version of the manuscript. The proportion of the funds invested in the j-th year of construction in the total investment r deb Rate of the debt fund r gov Discount rate of the public sector r pri Discount rate of the private sector e i The exterior benefits generated from project implementation in the i-th year Decision Variables α Initial investment for construction shared by the public sector. The private sector's share is notated as 1 − α. δ Subsidy per vehicle-kilometer ε Ratio of exterior benefits to be collected ρ Toll of per vehicle-kilometer τ Concession term In Table A1, o i and m i are expressed separately as follows.  NPV of the equity fund shared by the private sector S r gov NPV of the equity fund shared by the public sector F equity pri NPV of the subsidy paid by the public sector as costs, discounted by r gov irr pri Internal return rate (IRR) of the private sector over the project's life R sub NPV of the subsidy paid by the public sector and received by the private sector as revenue, discounted by r pri R tol NPV of toll collected by the private sector R ext NPV of the exterior benefits transferred and taken by the private sector as extra revenue C ope NPV of the operation costs paid by the private sector