Exemple d'utilisation de Python : Analyser les faiblesses d'un réseau¶
Ici nous allons utiliser les données suivantes : https://sergelhomme.fr/data/TD.py . Il faudra alors ouvrir Spyder pour ouvrir ce fichier dans l'éditeur et exécuter ce fichier comme ci-dessous. Le réseau s'affiche dans l'onglet "Graphes". Les différentes listes s'affichent dans l'onglet "Explorateur de variables". Désormais, sauf précision contraire, toutes les lignes de codes devront être écrites à la suite dans l'éditeur.
import math
import matplotlib.pyplot as plt
arcs = [(0,43), (0,56), (1,2), (1,6), (2,5), (2,14), (3,4), (3,17),(4,5), (4,16), (5,15), (5,14), (6,7), (6,14), (7,8), (7,14), (8,9), (8,11), (10,11), (10,24), (11,12), (11,24), (11,25), (12,13), (12,23), (13,14), (13,22), (14,21), (15,21), (15,18), (16,17), (16,18), (17,19), (18,20), (19,20), (20,21), (20,35), (20,36), (21,33), (22,26), (22,27), (23,25), (23,26), (24,25), (25,30), (26,31), (27,32), (27,33), (28,29), (28,30), (29,30), (29,43), (30,31), (32,42), (33,34), (35,38), (36,37), (36,38), (37,39), (38,39), (38,41), (38,40), (39,40), (39,47), (40,46), (34,41), (42,43), (42,44), (43,44), (44,55), (45,55), (45,52), (46,51), (47,48), (47,50), (48,49), (48,50), (49,50), (50,51), (51,52), (52,53), (53,54), (53,55), (55,56), (34,42), (34,45), (35,36)]
arcs_w = [(0,43,1), (0,56,1), (1,2,1), (1,6,1), (2,5,1), (2,14,1), (3,4,1), (3,17,1),(4,5,1), (4,16,1), (5,15,1), (5,14,1), (6,7,1), (6,14,1), (7,8,1), (7,14,1), (8,9,1), (8,11,1), (10,11,1), (10,24,1), (11,12,1), (11,24,1), (11,25,1), (12,13,1), (12,23,1), (13,14,1), (13,22,1), (14,21,1), (15,21,1), (15,18,1), (16,17,1), (16,18,1), (17,19,1), (18,20,1), (19,20,1), (20,21,1), (20,35,1), (20,36,1), (21,33,1), (22,26,1), (22,27,1), (23,25,1), (23,26,1), (24,25,1), (25,30,1), (26,31,1), (27,32,1), (27,33,1), (28,29,1), (28,30,1), (29,30,1), (29,43,1), (30,31,1), (32,42,1), (33,34,1), (35,38,1), (36,37,1), (36,38,1), (37,39,1), (38,39,1), (38,41,1), (38,40,1), (39,40,1), (39,47,1), (40,46,1), (34,41,1), (42,43,1), (42,44,1), (43,44,1), (44,55,1), (45,55,1), (45,52,1), (46,51,1), (47,48,1), (47,50,1), (48,49,1), (48,50,1), (49,50,1), (50,51,1), (51,52,1), (52,53,1), (53,54,1), (53,55,1), (55,56,1), (34,42,1), (34,45,1), (35,36,1)]
villes = ["Bayonne","Calais", "Lille", "Strasbourg", "Metz", "Reims", "Amiens", "Rouen", "Caen", "Cherbourg", "Brest", "Rennes", "Le Mans", "Chartres", "Paris", "Troyes", "Nancy", "Mulhouse", "Langres", "Besancon", "Dijon", "Sens", "Orleans", "Angers", "Vannes", "Nantes", "Tours", "Vierzon", "La Rochelle", "Saintes", "Niort", "Poitiers", "Limoges", "Vichy", "Clermont-Fd", "Macon", "Bourg-en-Bresse", "Geneve", "Lyon", "Grenoble", "Valence", "St-Etienne", "Brive", "Bordeaux", "Montauban", "Millau", "Avignon", "Digne-les-bains", "Nice", "Toulon", "Marseille", "Nimes", "Montpellier", "Narbonne", "Perpignan", "Toulouse", "Pau"]
xcoord = [290536, 566071, 651612, 999047, 880353 ,723240, 597143, 510403, 402458, 314252, 95121, 300873, 440320, 537130, 600217, 729191, 883209, 974902 ,824164, 878735, 804396, 670312, 567614, 382062, 217305, 305434, 475537, 579635, 330654, 368990, 384445, 446713, 516171, 684134, 658346, 792006, 822968, 883940, 795064, 866204, 801797, 760350, 536794, 369261, 521176, 659532, 798657, 912156, 997330, 892862, 846478, 762898, 724347, 654405, 645736, 527848, 379897]
ycoord = [1839901, 2661635, 2626710, 2412056, 2464886, 2475527, 2543973, 2494702, 2468337, 2522966, 2398724, 2353777, 2336045, 2383537, 2428986, 2368090, 2417370, 2317371, 2322392, 2255629, 2261516, 2355793, 2322615, 2278686, 2307396, 2253671, 2267165, 2246988, 2134800, 2087224, 2151008, 2177518, 2093169, 2125746, 2086868, 2148329, 2137972, 2144712, 2087442, 2026343, 1995763, 2050138, 2017905, 1986275, 1891190, 1899931, 1886158, 1906979, 1868053, 1798209, 1815584, 1872733, 1846609, 1798260, 1743848, 1845136, 1814694]
arcs_geom = []
arcs_dist = []
for i in range(len(arcs)):
arcs_geom = arcs_geom + [[[xcoord[arcs[i][0]],xcoord[arcs[i][1]]],[ycoord[arcs[i][0]],ycoord[arcs[i][1]]]]]
arcs_dist = arcs_dist + [(arcs[i][0], arcs[i][1], math.sqrt((xcoord[arcs[i][0]] - xcoord[arcs[i][1]]) * (xcoord[arcs[i][0]] - xcoord[arcs[i][1]]) + (ycoord[arcs[i][0]] - ycoord[arcs[i][1]]) * (ycoord[arcs[i][0]] - ycoord[arcs[i][1]])))]
for i in range(len(arcs_geom)):
plt.plot(arcs_geom[i][0],arcs_geom[i][1],c='gray')
plt.plot(xcoord,ycoord,'o')
plt.axis('equal')
plt.show()
Le graphique montre le réseau sur lequel on va travailler : un réseau simplifié du réseau autoroutier français. On dispose d'une liste d'arcs (arcs) qui permet de savoir les sommets qui sont reliés par un arc. Ainsi, la première liste de la liste arcs permet de savoir que le sommet 0 est relié au sommet 43. La liste villes permet d'utiliser les valeurs des sommets comme index pour connaitre les villes concernées. A partir du nom d'une ville on peut obtenir son index avec la fonction index(). La liste arcs_dist permet de connaitre la distance de ces relations.
La séquence ci-dessous peut être effectuée directement dans la console pour éviter d'alourdir le code.
print(arcs[0])
print(villes[0])
print(villes[43])
print(villes[arcs[0][0]]) #aternative plus subtile pour utiliser les valeurs de la liste arcs
print(villes[arcs[0][1]]) #aternative plus subtile pour utiliser les valeurs de la liste arcs
print(villes.index('Paris'))
print(arcs_dist[0])
print(str(arcs_dist[0][2]/1000) + ' km') #utilisation du + pour concatener du texte et afficher la distance entre Bayonne et Bordeaux
(0, 43) Bayonne Bordeaux Bayonne Bordeaux 14 (0, 43, 166201.6049892419) 166.2016049892419 km
Les débuts avec Networkx : création d'un graphe et calcul d'indicateurs locaux¶
L'import de la bibliothèque Networkx et la création d'un graphe G est très simple :
import networkx as nx
G = nx.Graph()
Si vous faites un print(G) vous verrez qu'un objet Graph() a bien été créé, mais celui-ci est vide. A noter que ce graphe est non dirigé (c'est à dire que l'on peut le "parcourir" dans les deux sens, par exemple la relation entre Bayonne et Bordeaux implique que l'on peut aller de Bordeaux à Bayonne et de Bordeaux à Bayonne). Il est possible de créer des Graphes dirigés avec DiGraph(). Pour ajouter des objets à votre graphes, vous pouvez utiliser les fonctions add_edges_from() pour les graphes non pondérés (comme par exemple la liste arcs) ou add_weighted_edges_from() pour les graphes pondérés (comme par exemple la liste arcs_dist qui contient les distances des arcs).
G.add_weighted_edges_from(arcs_dist)
Si vous faites un print(G), rien de nouveau. Néanmoins, si vous tapez dans la console G.edges() vous verrez la liste des arcs ajoutés. G.edges(data=True) permet de voir les pondérations (ici les distances). G.nodes() permet de récupérer la liste des noeuds. Vous pouvez alors calculer de nombreux indicateurs de graphes : https://networkx.org/documentation/stable/reference/algorithms/index.html. Ci-dessous, on calcule les degrés des sommets, leur centralité intermédiaire et leur proximité. Ce sont des indicateurs de centralité très utilisés pour identifier les sommets "principaux/centraux" d'un graphe.
deg = G.degree()
bet = nx.betweenness_centrality(G, weight='weight')
clos = nx.closeness_centrality(G, distance='weight')
Pour faire les calculs de manière topologique, sans tenir compte des distances, il suffit de retirer les arguments : weight='weight' ou distance='weight'. Le degré est un indicateur relativement pauvre, apportant peu d'informations, pour un réseau d'infrastructures. En revanche, centralité intermédiaire et proximité sont davantage liés.
import numpy as np
plt.plot(np.array(list(bet.items()))[:,1],np.array(list(clos.items()))[:,1],'o')
plt.xlabel('Betweenness')
plt.ylabel('Closeness')
plt.show()
Pour consulter les valeurs obtenues, vous pouvez aller dans l'onglet "Explorateur de variables", puis double-cliquer sur l'objet. Une table s'ouvre, vous pouvez trier les valeurs. La ville la plus centrale (centralité intermédiaire) est la 38, qui correspond à Lyon. Si on veut afficher les dix villes les plus centrales, on pourra trier par ordre décroissant la liste bet (qui au passage est un dictionnaire) et utiliser une boucle.
import operator
tribet = sorted(bet.items(),key=operator.itemgetter(1),reverse=True)
print('------Betweenness--------')
for i in range(10):
print( str(i+1) + ' : ' + villes[tribet[i][0]] )
------Betweenness-------- 1 : Lyon 2 : Paris 3 : Dijon 4 : Clermont-Fd 5 : Chartres 6 : Orleans 7 : Sens 8 : Macon 9 : St-Etienne 10 : Bordeaux
Calcul de plus courts chemins¶
A priori l'utilisation principale de Networkx qui lui vaut d'être installé par défaut dans la distribution d'Anaconda est que cette bibliothèque permet des calculs de plus courts chemins de flow minimal...
Cette partie "Calcul de plus courts chemins" peut être effectuée directement dans la console pour éviter d'alourdir le code.
La fonction shortest_path() permet de calculer un plus court chemin entre deux sommets et retourne les sommets empruntés par ce plus court chemin. La fonction shortest_path_length() fonctionne de manière identique mais retourne la longueur de ce plus court chemin. L'argument weight='weight' permet de tenir compte des poids/distances des arcs.
print(nx.shortest_path(G,10,12))
print(nx.shortest_path_length(G,10,12))
print(nx.shortest_path(G,10,12,weight='weight'))
print(nx.shortest_path_length(G,10,12,weight='weight'))
[10, 11, 12] 2 [10, 11, 12] 351174.0490480197
Pour calculer la longueur de tous les plus court chemins, shortest_path_length() peut être utilisé sans préciser de sommets et produit un dictionnaire que l'on peut parcourir.
spl = dict(nx.shortest_path_length(G,weight='weight'))
print(spl[10][12])
print(spl[34][0])
351174.0490480197 476446.89317854075
Pour calculer la somme des distances d'une ville (son éloignement) ou l'inverse des distances (sa proximité), on peut effectuer une boucle. Pour l'inverse de la distance et éviter une division par zéro il faudra rajouter une condition.
el = 0
for i in range(57) :
el = el + spl[34][i]
print(el)
22309106.083352223
prox = 0
for i in range(57) :
if i != 14 :
prox = prox + (1/spl[14][i])
print(prox)
0.00018418207796682933
Pour obtenir ces valeurs sur l'ensemble du graphe, il faudra alors effectuer une double boucle qui parcourt tous les sommets et pas seulement Paris (sommet 14) comme ci-dessus. Il est aussi possible d'utiliser la fonction average_shortest_path_length(G,weight='weight') qui moyenne cette valeur pour son éloignement. Malheureusement, la même fonction n'existe pas pour la proximité.
aspl = nx.average_shortest_path_length(G,weight='weight')
print(aspl)
516309.01461460255
Cette valeur moyenne peut être interprétée comme un indicateur de performance du réseau, plus exactement comme un indicateur. En effet, plus le réseau est connecté par des routes "courtes" pour relier les villes, plus cette valeur moyenne baissera et le réseau sera efficace (moins de distances parcourues par les usagers, moins de pollutions, moins de coûts de transports...). Dans le cas contraire, cette moyenne augmentera et le réseau sera moins efficace. On peut alors supprimer un sommet à l'aide de la fonction remove_node() et étudier les conséquences. Ici un nouveau graphe est créé (Gper) pour éviter de toucher au graphe initial (non perturbé).
Gper = nx.Graph()
Gper.add_weighted_edges_from(arcs_dist)
Gper.remove_node(14)
asplper = nx.average_shortest_path_length(Gper,weight='weight')
print(asplper - aspl)
59427.51575177431
En moyenne, si vous supprimez Paris (sommet 14) du graphe, vous augmentez les plus court chemins (les parcours) de 59 km.
Etudier la fragilité du réseau¶
Pour étudier la fragilité du réseau, on peut alors simuler toutes les suppressions de tous les sommets à l'aide d'une boucle, en créant à chaque itération un nouveau graphe à perturber (G2). On stocke les résultats (identifiant du sommet et valeur) dans une liste nommée vul.
aspl = nx.average_shortest_path_length(G,weight='weight')
vul = []
for i in range(57):
if i != 8 and i!=53 :
G2 = nx.Graph()
G2.add_weighted_edges_from(arcs_dist)
G2.remove_node(i)
aspl2 = nx.average_shortest_path_length(G2,weight='weight')
v = aspl2 - aspl
vul = vul + [[i,v]]
La suppression des sommets 8 (Caen) et 53 (Narbonne) pose problème, car il coupe le graphe en différentes composantes connexes en isolant respectivement Cherbourg et Perpignant. C'est pour cela qu'il faut créer des exceptions. L'explorateur de variables permet de voir quels sont les sommets qui causent le plus de détours. Néanmoins, on peut trier la liste vul et afficher les 10 villes causant le plus de détours dans la console.
trivul = sorted(vul,key=operator.itemgetter(1),reverse=True)
print('------Vulnerabilite--------')
for i in range(10):
print( str(i+1) + ' : ' + villes[trivul[i][0]] )
------Vulnerabilite-------- 1 : Paris 2 : Dijon 3 : Clermont-Fd 4 : Chartres 5 : Lyon 6 : Niort 7 : Bordeaux 8 : Montpellier 9 : Orleans 10 : Saintes
Pour tenir compte des sommets 8 (Caen) et 53 (Narbonne) et ainsi tenir compte des pertes de connectivité induites par la suppression de ces sommets, il faut effectuer les calculs en termes de proximité (inverse des distances) et non en termes d'éloignement. Cela évite les valeurs infinies. On peut alors créer notre propre fonction eff_total() qui effectue ce calcul.
def eff_total(G):
spl = dict(nx.shortest_path_length(G,weight='weight'))
sspl = 0
for i in range(57):
for j in range(57):
if i != j :
try :
sspl = sspl + 1 / spl[i][j]
except :
sspl = sspl + 0
return sspl
sspl = eff_total(G)
print(sspl)
0.008854532598888683
La fonction ci-dessus utilise les tests, try : except : , dans le but d'éviter les erreurs causées par l'absence de valeurs de distance calculées lors des pertes de connectivité. On peut alors effectuer un boucle, en créant à chaque itération un nouveau graphe à perturber (G2), que va calculer les perturbations engendrées dans un objet nommé ici vul2.
vul2 =[]
for k in range(57):
G2 = nx.Graph()
G2.add_weighted_edges_from(arcs_dist)
G2.remove_node(k)
sspl2 = eff_total(G2)
vul2 = vul2 + [[k,sspl - sspl2]]
Comme précédemment, on peut trier la liste vul2 et afficher les 10 villes causant le plus d'impacts dans la console.
trivul2 = sorted(vul2,key=operator.itemgetter(1),reverse=True)
print('------Vulnerabilite2--------')
for i in range(10):
print( str(i+1) + ' : ' + villes[trivul2[i][0]] )
------Vulnerabilite2-------- 1 : Paris 2 : Lyon 3 : Dijon 4 : Nimes 5 : Chartres 6 : Clermont-Fd 7 : Montpellier 8 : Niort 9 : Narbonne 10 : Orleans
Ce deuxième indicateur modifie le premier classement établi et permet surtout d'avoir une valeur de vulnérabilité pour chaque sommet. Ainsi, il est possible de tester des corrélations entre la vulnérabilité et les valeurs de centralité.
tribet = sorted(list(bet.items()))
atribet = np.array(tribet)
avul2 = np.array(vul2)
plt.plot(atribet[:,1],avul2[:,1],'o')
plt.xlabel('Betweenness')
plt.ylabel('Vulnerability')
plt.show()
import scipy.stats as sts
reg = sts.linregress(atribet[:,1],avul2[:,1])
print("Le r2 est de : " + str(reg[2]*reg[2]))
Le r2 est de : 0.7926181245406703
Bonus : dans l'intervalle de confiance ou pas ?¶
yapprox = reg[0] * atribet[:,1] + reg[1]
t = sts.t.ppf(0.975,len(G.nodes())-2) #risque0.05
moy = np.mean(atribet[:,1])
v = np.var(atribet[:,1])
n = len(G.nodes())
su = reg[4] * math.sqrt(n*v)
ytolere = []
ytolere2 = []
for i in range(n):
val = yapprox[i] + (t*su*math.sqrt(1+(1/float(n))+((atribet[i,1]-moy)*(atribet[i,1]-moy)/float(n*v))))
ytolere = ytolere + [val]
val2 = yapprox[i] - (t*su*math.sqrt(1+(1/float(n))+((atribet[i,1]-moy)*(atribet[i,1]-moy)/float(n*v))))
ytolere2 = ytolere2 + [val2]
#plt.subplots(figsize=(6, 4))
fig = plt.figure()
plt.plot(atribet[:,1],avul2[:,1],'o', markersize = 3)
plt.plot(sorted(atribet[:,1]),sorted(yapprox), color = 'r',linewidth=1)
plt.plot(sorted(atribet[:,1]),sorted(ytolere), '--',color = 'orange',linewidth=1)
plt.plot(sorted(atribet[:,1]),sorted(ytolere2), '--', color = 'orange',linewidth=1)
plt.xlabel('Betweenness')
plt.ylabel('Vulnerability')
plt.show()
Les deux points les plus éloignés de l'intervalle de confiance sont Caen et Narbonne. Si le passage par l'inverse des distances (proximité) permet ici de mélanger les problématiques de pertes d'efficacité (éloignement des trajets) et pertes de connectivité (impossibilité de faire le trajet), cela relève plus de l'astuce mathématique. Pour faire rentrer dans le rang ces deux sommets on peut jouer sur la valeur de la puissance à utiliser...