|
8 | 8 | % With all_routes, need to set beta = 1 and CrossGoodCongestion 'off' to use efficient dual solution |
9 | 9 |
|
10 | 10 | % Read Undirected Graph |
11 | | -graph = readtable("data/transport_network/trans_african/fastest_routes_graph_edges.csv"); |
12 | | -% histogram(graph.ug_cost, 'BinWidth', 100); |
| 11 | +edges = readtable("data/transport_network/trans_african/fastest_routes_graph_edges.csv"); |
| 12 | +% histogram(edges.ug_cost, 'BinWidth', 100); |
13 | 13 | % Adjusting |
14 | | -graph.distance = graph.distance / 1000; % Convert to km |
15 | | -graph.border_dist = graph.border_dist / 1000; % Convert to km |
16 | | -graph.duration = graph.duration / 60; % Convert to hours |
17 | | -graph.total_cost = graph.ug_cost / 1e6; % Convert to millions |
| 14 | +edges.distance = edges.distance / 1000; % Convert to km |
| 15 | +edges.border_dist = edges.border_dist / 1000; % Convert to km |
| 16 | +edges.duration = edges.duration / 60; % Convert to hours |
| 17 | +edges.total_cost = edges.ug_cost / 1e6; % Convert to millions |
18 | 18 |
|
19 | | -n = max(max(graph.from), max(graph.to)); |
| 19 | +n = max(max(edges.from), max(edges.to)); |
20 | 20 |
|
21 | 21 | % Create Adjacency Matrix |
22 | 22 | adj_matrix = false(n, n); |
23 | | -for i = 1:height(graph) |
24 | | - adj_matrix(graph.from(i), graph.to(i)) = true; |
25 | | - adj_matrix(graph.to(i), graph.from(i)) = true; |
| 23 | +for i = 1:height(edges) |
| 24 | + adj_matrix(edges.from(i), edges.to(i)) = true; |
| 25 | + adj_matrix(edges.to(i), edges.from(i)) = true; |
26 | 26 | end |
27 | 27 |
|
28 | 28 | % Read Nodes Data |
|
33 | 33 |
|
34 | 34 | % Create Infrastructure Matrix: Following Graff (2024) = average speed in km/h: length of route is accounted for in cost function |
35 | 35 | infra_matrix = zeros(n, n); |
36 | | -for i = 1:height(graph) |
37 | | - speed = graph.distance(i) / graph.duration(i); |
38 | | - infra_matrix(graph.from(i), graph.to(i)) = speed; |
39 | | - infra_matrix(graph.to(i), graph.from(i)) = speed; |
| 36 | +for i = 1:height(edges) |
| 37 | + speed = edges.distance(i) / edges.duration(i); |
| 38 | + infra_matrix(edges.from(i), edges.to(i)) = speed; |
| 39 | + infra_matrix(edges.to(i), edges.from(i)) = speed; |
40 | 40 | end |
41 | | -summary(graph.distance ./ graph.duration); |
| 41 | +summary(edges.distance ./ edges.duration); |
42 | 42 |
|
43 | 43 | % Create Iceberg Trade Cost Matrix: Following Graff (2024) |
44 | | -% graph.distance = graph.distance + graph.border_dist; % uncomment to enable frictions |
| 44 | +% edges.distance = edges.distance + edges.border_dist; % uncomment to enable frictions |
45 | 45 | iceberg_matrix = zeros(n, n); |
46 | | -for i = 1:height(graph) |
47 | | - iceberg_matrix(graph.from(i), graph.to(i)) = 0.1158826 * log(graph.distance(i) / 1.609); |
48 | | - iceberg_matrix(graph.to(i), graph.from(i)) = 0.1158826 * log(graph.distance(i) / 1.609); |
| 46 | +for i = 1:height(edges) |
| 47 | + iceberg_matrix(edges.from(i), edges.to(i)) = 0.1158826 * log(edges.distance(i) / 1.609); |
| 48 | + iceberg_matrix(edges.to(i), edges.from(i)) = 0.1158826 * log(edges.distance(i) / 1.609); |
49 | 49 | end |
50 | 50 | iceberg_matrix(iceberg_matrix < 0) = 0; |
51 | 51 |
|
52 | 52 | % Create Infrastructure Building Cost Matrix: Following Graff (2024) |
53 | 53 | infra_building_matrix = zeros(n, n); |
54 | | -for i = 1:height(graph) |
55 | | - infra_building_matrix(graph.from(i), graph.to(i)) = graph.total_cost(i); |
56 | | - infra_building_matrix(graph.to(i), graph.from(i)) = graph.total_cost(i); |
| 54 | +for i = 1:height(edges) |
| 55 | + infra_building_matrix(edges.from(i), edges.to(i)) = edges.total_cost(i); |
| 56 | + infra_building_matrix(edges.to(i), edges.from(i)) = edges.total_cost(i); |
57 | 57 | end |
58 | 58 |
|
59 | 59 | % Basic characteristics of the economy |
|
110 | 110 | 'a', a, 'sigma', sigma, 'N', N, 'alpha', alpha, 'beta', beta, 'gamma', gamma, 'rho', rho, ... |
111 | 111 | 'verbose', 'on', 'K', K, 'TolKappa', 1e-5); |
112 | 112 |
|
113 | | -[param, g] = create_graph(param,[],[],'X',nodes.lon,'Y',nodes.lat,'Type','custom','Adjacency',adj_matrix); |
| 113 | +[param, graph] = create_graph(param,[],[],'X',nodes.lon,'Y',nodes.lat,'Type','custom','Adjacency',adj_matrix); |
114 | 114 |
|
115 | 115 | param.Lj = population; |
116 | 116 | param.Zjn = productivity; |
117 | 117 | % param.omegaj = weights; % In case want to give more weight to specific locations |
118 | 118 | param.Hj = population .* (1-alpha); % TG: I normalise this because the general utility function has a (h_j/(1-alpha))^(1-alpha) thing with it |
119 | 119 |
|
120 | | -g.delta_i = infra_building_matrix; |
121 | | -g.delta_tau = iceberg_matrix; |
| 120 | +graph.delta_i = infra_building_matrix; |
| 121 | +graph.delta_tau = iceberg_matrix; |
122 | 122 |
|
123 | 123 | % Naming conventions: if IRS -> 'cgc_irs' or 'irs_na' without annealing; if alpha = 0.1 -> add '_alpha01'; if with_ports = false -> add '_noport'; if frictions -> add '_bc' (border cost) |
124 | 124 | filename = '22g_10b_fixed_cgc_sigma38_rho0'; % adjust if sigma != 1.5 |
125 | 125 | fprintf('File extension: %s\n', filename); |
126 | 126 |
|
127 | 127 | % Solve allocation from existing infrastructure |
128 | 128 | strcat("Started P_stat on ", datestr(datetime('now'))) |
129 | | -res_stat = solve_allocation(param, g, infra_matrix, true); |
| 129 | +res_stat = solve_allocation(param, graph, infra_matrix, true); |
130 | 130 |
|
131 | 131 | % Solve Optimal Network (this can take long - up to 48h) |
132 | 132 | strcat("Started P_opt on ", datestr(datetime('now'))) |
133 | | -res_opt = optimal_network(param, g, infra_matrix, min_mask, max_mask, false); |
| 133 | +res_opt = optimal_network(param, graph, infra_matrix, min_mask, max_mask, false); |
134 | 134 |
|
135 | 135 | % Check: should be 1 |
136 | | -K_ratio = K / sum(sum(res_opt.Ijk .* g.delta_i)); |
| 136 | +K_ratio = K / sum(sum(res_opt.Ijk .* graph.delta_i)); |
137 | 137 | disp(K_ratio); |
138 | 138 |
|
139 | 139 | % Plot Network |
140 | | -plot_graph(param, g, res_opt.Ijk, 'Margin', 0.01); |
| 140 | +plot_graph(param, graph, res_opt.Ijk, 'Margin', 0.01); |
141 | 141 | axis equal; |
142 | | -plot_graph(param, g, res_opt.Ijk - infra_matrix, 'Margin', 0.01); |
| 142 | +plot_graph(param, graph, res_opt.Ijk - infra_matrix, 'Margin', 0.01); |
143 | 143 | axis equal; |
144 | 144 |
|
145 | 145 | % Saving: Nodes |
|
164 | 164 | writetable(res_nodes, sprintf('results/transport_network/GE/trans_african/nodes_results_%s.csv', filename)); |
165 | 165 |
|
166 | 166 | % Saving: Graph / Edges |
167 | | -res_graph = graph; |
168 | | -res_graph.Ijk_orig = res_to_vec(infra_matrix, graph); |
169 | | -res_graph.Ijk = res_to_vec(res_opt.Ijk, graph); |
| 167 | +res_edges = edges; |
| 168 | +res_edges.Ijk_orig = res_to_vec(infra_matrix, edges); |
| 169 | +res_edges.Ijk = res_to_vec(res_opt.Ijk, edges); |
170 | 170 | for n = 1:N |
171 | | - res_graph = setfield(res_graph, ['Qjk_', num2str(n)], res_to_vec(res_opt.Qjkn(:,:,n), graph)); |
| 171 | + res_edges = setfield(res_edges, ['Qjk_', num2str(n)], res_to_vec(res_opt.Qjkn(:,:,n), edges)); |
172 | 172 | end |
173 | | -writetable(res_graph, sprintf('results/transport_network/GE/trans_african/edges_results_%s.csv', filename)); |
| 173 | +writetable(res_edges, sprintf('results/transport_network/GE/trans_african/edges_results_%s.csv', filename)); |
0 commit comments