- Implement Dijkstra's shortest path algorithm in GAMS,
- implement Dijkstra's shortest path algorithm in Python and
- use a simple LP model.
Data
*-----------------------------------------------------* set up network*-----------------------------------------------------set i 'nodes' /node1*node50/;alias (i,j,v);set a(i,j) 'arcs of sparse network';* no self loops and density of 15%a(i,j)$(ord(i)<>ord(j)) = uniform(0,1)<0.15;singleton setssource(i) /node1/sink(i) /node49/;* sink is chosen so optimal path needs quite a few hopsparameter cost(i,j) 'random costs (lengths)';cost(a) = uniform(1,10);option a:0:0:7, cost:3:0:6;display$(card(i)<=50)"============== data ==============",i,source,sink,a,cost;
- Singleton sets can be used to enforce that the set only contains zero or one element.
- The option statements produce a more compact display output.
- The $ condition on the display prevents printing for very large data sets. It is better to use the GDX data viewer for those.
Listing
---- 36 ============== data ============== ---- 36 SET i nodes node1 , node2 , node3 , node4 , node5 , node6 , node7 , node8 , node9 , node10, node11, node12 node13, node14, node15, node16, node17, node18, node19, node20, node21, node22, node23, node24 node25, node26, node27, node28, node29, node30, node31, node32, node33, node34, node35, node36 node37, node38, node39, node40, node41, node42, node43, node44, node45, node46, node47, node48 node49, node50 ---- 36 SET source node1 ---- 36 SET sink node49 ---- 36 SET a arcs of sparse network node1 .node10, node1 .node16, node1 .node24, node1 .node32, node1 .node43, node1 .node45, node2 .node7 node2 .node8 , node2 .node11, node2 .node13, node2 .node18, node2 .node25, node2 .node30, node2 .node32 node2 .node34, node2 .node49, node3 .node1 , node3 .node4 , node3 .node5 , node3 .node15, node3 .node22 node3 .node24, node3 .node28, node3 .node31, node3 .node43, node4 .node6 , node4 .node18, node4 .node20 node4 .node22, node4 .node24, node4 .node29, node4 .node35, node4 .node37, node4 .node38, node4 .node41 node5 .node8 , node5 .node14, node5 .node22, node5 .node26, node5 .node32, node5 .node38, node5 .node48 node5 .node50, node6 .node9 , node6 .node16, node6 .node20, node6 .node23, node6 .node27, node6 .node34 node6 .node39, node6 .node43, node6 .node46, node6 .node50, node7 .node11, node7 .node26, node7 .node38 node7 .node39, node7 .node44, node7 .node45, node8 .node1 , node8 .node4 , node8 .node15, node8 .node18 node8 .node20, node8 .node28, node8 .node33, node8 .node46, node8 .node49, node9 .node1 , node9 .node11 node9 .node26, node9 .node50, node10.node2 , node10.node17, node10.node27, node10.node36, node10.node38 node10.node41, node10.node48, node10.node50, node11.node14, node11.node31, node11.node34, node11.node50 node12.node2 , node12.node4 , node12.node9 , node12.node10, node12.node17, node12.node39, node12.node42 node12.node43, node13.node4 , node13.node7 , node13.node12, node13.node23, node13.node24, node13.node33 node13.node41, node13.node45, node13.node49, node13.node50, node14.node2 , node14.node3 , node14.node5 node14.node15, node14.node18, node14.node21, node14.node39, node14.node45, node14.node47, node14.node50 node15.node10, node15.node13, node15.node16, node15.node22, node15.node31, node15.node38, node16.node1 node16.node6 , node16.node8 , node16.node10, node16.node44, node17.node4 , node17.node6 , node17.node13 node17.node19, node17.node22, node17.node24, node17.node35, node17.node36, node17.node37, node18.node4 node18.node5 , node18.node12, node18.node13, node18.node37, node18.node41, node18.node42, node18.node45 node19.node1 , node19.node10, node19.node14, node19.node18, node19.node31, node19.node34, node19.node50 node20.node9 , node20.node11, node20.node16, node20.node19, node20.node22, node20.node23, node20.node26 node20.node27, node20.node29, node20.node30, node20.node41, node20.node46, node20.node49, node21.node1 node21.node5 , node21.node11, node21.node23, node21.node29, node21.node31, node21.node46, node22.node9 node22.node12, node22.node18, node22.node21, node22.node29, node22.node39, node22.node50, node23.node3 node23.node8 , node23.node17, node23.node19, node23.node21, node23.node35, node23.node37, node23.node41 node23.node50, node24.node7 , node24.node9 , node24.node10, node24.node18, node24.node44, node24.node45 node24.node47, node25.node5 , node25.node6 , node25.node32, node25.node34, node25.node45, node25.node47 node26.node3 , node26.node11, node26.node13, node26.node14, node26.node15, node26.node19, node26.node28 node26.node33, node26.node36, node26.node39, node26.node41, node26.node44, node26.node48, node27.node8 node27.node10, node27.node12, node27.node15, node27.node49, node28.node9 , node28.node10, node28.node13 node28.node15, node28.node27, node28.node47, node29.node13, node29.node21, node29.node24, node29.node46 node30.node2 , node30.node6 , node30.node16, node30.node25, node30.node32, node30.node33, node30.node34 node30.node36, node30.node42, node30.node46, node30.node48, node31.node16, node31.node18, node31.node19 node31.node20, node31.node29, node31.node32, node31.node33, node32.node2 , node32.node4 , node32.node8 node32.node9 , node32.node11, node32.node18, node32.node22, node32.node29, node32.node31, node32.node38 node32.node45, node32.node47, node32.node49, node33.node8 , node33.node11, node33.node12, node33.node22 node33.node29, node33.node30, node33.node31, node33.node41, node33.node48, node34.node20, node34.node22 node34.node47, node35.node1 , node35.node4 , node35.node14, node35.node21, node35.node23, node35.node28 node35.node29, node35.node30, node35.node31, node35.node40, node35.node43, node36.node4 , node36.node10 node36.node13, node36.node20, node36.node26, node36.node34, node36.node43, node37.node2 , node37.node4 node37.node8 , node37.node12, node37.node13, node37.node19, node37.node22, node37.node30, node37.node39 node37.node40, node37.node49, node38.node2 , node38.node3 , node38.node4 , node38.node7 , node38.node12 node38.node13, node38.node17, node38.node44, node38.node49, node39.node6 , node39.node11, node39.node17 node39.node27, node39.node34, node39.node45, node39.node49, node40.node5 , node40.node15, node40.node21 node40.node24, node40.node26, node40.node38, node40.node50, node41.node18, node41.node23, node41.node33 node41.node38, node41.node43, node41.node46, node42.node1 , node42.node14, node42.node34, node42.node38 node43.node3 , node43.node6 , node43.node16, node43.node24, node43.node25, node43.node34, node43.node37 node43.node40, node43.node48, node44.node4 , node44.node8 , node44.node11, node44.node20, node44.node22 node44.node24, node44.node34, node44.node43, node44.node46, node45.node4 , node45.node10, node45.node18 node45.node28, node45.node29, node45.node35, node45.node41, node45.node44, node45.node46, node46.node2 node46.node6 , node46.node7 , node46.node14, node46.node15, node46.node21, node46.node22, node46.node26 node46.node29, node46.node33, node46.node41, node46.node47, node47.node7 , node47.node12, node47.node19 node47.node28, node47.node31, node47.node36, node47.node39, node48.node3 , node48.node5 , node48.node8 node48.node9 , node48.node29, node49.node6 , node49.node21, node49.node22, node49.node41, node49.node47 node50.node8 , node50.node9 , node50.node13, node50.node19, node50.node24, node50.node29, node50.node34 node50.node39, node50.node41, node50.node47 ---- 36 PARAMETER cost random costs (lengths) node1 .node10 4.935, node1 .node16 4.958, node1 .node24 7.343, node1 .node32 3.319, node1 .node43 6.534 node1 .node45 8.792, node2 .node7 2.969, node2 .node8 4.780, node2 .node11 3.507, node2 .node13 8.732 node2 .node18 4.993, node2 .node25 5.997, node2 .node30 1.826, node2 .node32 2.595, node2 .node34 5.995 node2 .node49 3.239, node3 .node1 6.640, node3 .node4 8.003, node3 .node5 3.465, node3 .node15 8.997 node3 .node22 6.381, node3 .node24 1.738, node3 .node28 2.199, node3 .node31 9.425, node3 .node43 8.258 node4 .node6 6.076, node4 .node18 2.368, node4 .node20 7.365, node4 .node22 9.519, node4 .node24 5.974 node4 .node29 2.300, node4 .node35 7.104, node4 .node37 2.546, node4 .node38 9.192, node4 .node41 4.910 node5 .node8 8.590, node5 .node14 8.442, node5 .node22 6.276, node5 .node26 6.035, node5 .node32 5.902 node5 .node38 9.634, node5 .node48 9.818, node5 .node50 3.369, node6 .node9 1.503, node6 .node16 6.245 node6 .node20 5.675, node6 .node23 1.090, node6 .node27 6.286, node6 .node34 1.174, node6 .node39 4.131 node6 .node43 2.046, node6 .node46 8.588, node6 .node50 8.779, node7 .node11 5.909, node7 .node26 1.352 node7 .node38 2.318, node7 .node39 5.750, node7 .node44 4.929, node7 .node45 7.365, node8 .node1 8.236 node8 .node4 7.909, node8 .node15 1.318, node8 .node18 6.358, node8 .node20 6.210, node8 .node28 5.611 node8 .node33 3.718, node8 .node46 9.171, node8 .node49 5.474, node9 .node1 6.001, node9 .node11 4.602 node9 .node26 2.126, node9 .node50 2.395, node10.node2 5.468, node10.node17 1.925, node10.node27 4.020 node10.node36 5.135, node10.node38 2.461, node10.node41 9.983, node10.node48 2.744, node10.node50 3.616 node11.node14 1.502, node11.node31 3.839, node11.node34 7.667, node11.node50 7.526, node12.node2 9.957 node12.node4 7.407, node12.node9 8.228, node12.node10 2.872, node12.node17 6.631, node12.node39 4.910 node12.node42 9.873, node12.node43 6.506, node13.node4 7.257, node13.node7 6.780, node13.node12 2.237 node13.node23 4.215, node13.node24 2.010, node13.node33 2.788, node13.node41 4.470, node13.node45 6.519 node13.node49 7.108, node13.node50 4.866, node14.node2 4.525, node14.node3 5.524, node14.node5 5.575 node14.node15 8.529, node14.node18 9.698, node14.node21 2.993, node14.node39 6.016, node14.node45 4.893 node14.node47 8.876, node14.node50 9.674, node15.node10 3.561, node15.node13 10.000, node15.node16 8.717 node15.node22 2.936, node15.node31 9.018, node15.node38 9.255, node16.node1 4.665, node16.node6 6.419 node16.node8 8.808, node16.node10 1.581, node16.node44 7.033, node17.node4 9.726, node17.node6 8.567 node17.node13 6.224, node17.node19 2.761, node17.node22 5.650, node17.node24 8.556, node17.node35 3.517 node17.node36 8.890, node17.node37 6.786, node18.node4 9.892, node18.node5 6.459, node18.node12 3.442 node18.node13 2.695, node18.node37 4.357, node18.node41 1.340, node18.node42 5.210, node18.node45 9.967 node19.node1 5.447, node19.node10 7.270, node19.node14 3.699, node19.node18 6.909, node19.node31 2.578 node19.node34 8.683, node19.node50 4.528, node20.node9 5.646, node20.node11 5.753, node20.node16 9.444 node20.node19 4.427, node20.node22 7.112, node20.node23 2.080, node20.node26 3.944, node20.node27 4.051 node20.node29 5.321, node20.node30 9.880, node20.node41 5.980, node20.node46 1.884, node20.node49 8.807 node21.node1 9.482, node21.node5 7.327, node21.node11 3.316, node21.node23 9.991, node21.node29 2.960 node21.node31 1.958, node21.node46 9.173, node22.node9 9.722, node22.node12 3.332, node22.node18 2.191 node22.node21 2.978, node22.node29 2.612, node22.node39 2.027, node22.node50 2.793, node23.node3 9.788 node23.node8 4.898, node23.node17 1.708, node23.node19 1.526, node23.node21 3.468, node23.node35 2.787 node23.node37 4.987, node23.node41 8.546, node23.node50 4.640, node24.node7 6.768, node24.node9 9.047 node24.node10 3.371, node24.node18 9.678, node24.node44 6.026, node24.node45 2.934, node24.node47 7.918 node25.node5 8.576, node25.node6 9.978, node25.node32 3.333, node25.node34 7.252, node25.node45 6.123 node25.node47 9.441, node26.node3 3.345, node26.node11 6.221, node26.node13 1.543, node26.node14 9.585 node26.node15 2.702, node26.node19 2.517, node26.node28 8.481, node26.node33 1.996, node26.node36 1.350 node26.node39 1.547, node26.node41 4.792, node26.node44 9.055, node26.node48 6.946, node27.node8 6.228 node27.node10 8.287, node27.node12 5.337, node27.node15 1.861, node27.node49 4.357, node28.node9 4.119 node28.node10 4.379, node28.node13 8.007, node28.node15 7.958, node28.node27 4.844, node28.node47 6.006 node29.node13 4.925, node29.node21 7.478, node29.node24 6.365, node29.node46 2.992, node30.node2 6.980 node30.node6 3.011, node30.node16 7.165, node30.node25 4.757, node30.node32 9.094, node30.node33 3.880 node30.node34 3.320, node30.node36 3.025, node30.node42 6.998, node30.node46 7.823, node30.node48 7.034 node31.node16 5.542, node31.node18 5.119, node31.node19 7.528, node31.node20 5.688, node31.node29 6.103 node31.node32 9.320, node31.node33 5.357, node32.node2 6.146, node32.node4 6.954, node32.node8 7.033 node32.node9 2.089, node32.node11 5.756, node32.node18 8.419, node32.node22 6.365, node32.node29 3.088 node32.node31 9.529, node32.node38 2.097, node32.node45 7.319, node32.node47 5.156, node32.node49 9.927 node33.node8 2.974, node33.node11 5.628, node33.node12 1.758, node33.node22 9.584, node33.node29 7.477 node33.node30 8.211, node33.node31 2.058, node33.node41 8.733, node33.node48 2.464, node34.node20 6.887 node34.node22 3.249, node34.node47 1.145, node35.node1 7.289, node35.node4 1.621, node35.node14 7.487 node35.node21 5.222, node35.node23 4.914, node35.node28 1.888, node35.node29 4.035, node35.node30 9.912 node35.node31 2.248, node35.node40 8.902, node35.node43 6.387, node36.node4 9.659, node36.node10 9.728 node36.node13 6.889, node36.node20 7.723, node36.node26 2.412, node36.node34 3.873, node36.node43 3.336 node37.node2 3.735, node37.node4 1.323, node37.node8 5.482, node37.node12 9.341, node37.node13 6.038 node37.node19 5.192, node37.node22 4.572, node37.node30 3.735, node37.node39 6.632, node37.node40 9.157 node37.node49 4.092, node38.node2 7.071, node38.node3 1.809, node38.node4 6.093, node38.node7 3.677 node38.node12 2.503, node38.node13 8.227, node38.node17 1.710, node38.node44 9.299, node38.node49 9.717 node39.node6 8.422, node39.node11 2.905, node39.node17 9.507, node39.node27 4.525, node39.node34 8.687 node39.node45 8.611, node39.node49 1.827, node40.node5 3.564, node40.node15 9.269, node40.node21 1.924 node40.node24 1.031, node40.node26 8.250, node40.node38 1.384, node40.node50 2.610, node41.node18 8.303 node41.node23 6.795, node41.node33 3.430, node41.node38 1.928, node41.node43 8.944, node41.node46 9.587 node42.node1 7.529, node42.node14 7.709, node42.node34 7.065, node42.node38 3.371, node43.node3 2.107 node43.node6 9.192, node43.node16 7.813, node43.node24 1.631, node43.node25 2.240, node43.node34 4.364 node43.node37 4.801, node43.node40 8.194, node43.node48 1.315, node44.node4 6.554, node44.node8 8.180 node44.node11 3.808, node44.node20 8.348, node44.node22 9.854, node44.node24 8.583, node44.node34 6.454 node44.node43 7.715, node44.node46 6.361, node45.node4 2.677, node45.node10 7.900, node45.node18 7.693 node45.node28 7.206, node45.node29 4.822, node45.node35 1.664, node45.node41 8.805, node45.node44 3.055 node45.node46 9.431, node46.node2 6.782, node46.node6 2.903, node46.node7 1.261, node46.node14 5.580 node46.node15 5.191, node46.node21 6.824, node46.node22 4.197, node46.node26 2.120, node46.node29 8.157 node46.node33 4.546, node46.node41 6.742, node46.node47 3.792, node47.node7 1.859, node47.node12 4.877 node47.node19 6.375, node47.node28 4.129, node47.node31 1.647, node47.node36 3.583, node47.node39 3.284 node48.node3 9.610, node48.node5 6.638, node48.node8 5.598, node48.node9 5.516, node48.node29 3.071 node49.node6 2.446, node49.node21 8.807, node49.node22 9.038, node49.node41 5.136, node49.node47 4.231 node50.node8 6.081, node50.node9 6.310, node50.node13 3.404, node50.node19 1.530, node50.node24 5.064 node50.node29 5.085, node50.node34 9.772, node50.node39 8.410, node50.node41 7.047, node50.node47 4.656
Dijkstra's algorithm
The method was actually already implemented years before [2]:
"Dijkstra formulated and solved the shortest path problem for a demonstration at the official inauguration of the ARMAC computer in 1956. Because of the absence of journals dedicated to automatic computing, he did not publish the result until 1959".
This returns two arrays. dist[i] is the optimal total distance (cost) from the source to node \(i\). prev[i] is the previous node on the optimal shortest path. To form a path \(S\) from this prev array, the following algorithm is proposed:
This is not the most sophisticated implementation. Real fast versions use data structures like a priority queue. But for our purposes, simplicity is more important than raw speed.
Approach 1: GAMS Implementation of Dijkstra
*-----------------------------------------------------* data structures*-----------------------------------------------------setsprev(i,j) 'i is previous wrt j'Q(i) 'unvisited nodes';parameter dist(i) 'distance from source';* initializeQ(i) = yes;dist(i) = INF;dist(source) = 0;*-----------------------------------------------------* algorithm*-----------------------------------------------------scalar mindist;singleton set u(i);scalar alt;* if multiple assignments to singleton set,* only the last one will remain. This will* guarantee a singleton set will hold* 0 or 1 elements.option strictSingleton=0;while(card(Q)>0,* find node in Q with minimum distance from the sourcemindist = smin(q,dist(q));u(Q) = dist(Q)=mindist;* remove u from QQ(u) = no;* loop over neighbors of u in Qloop(a(u,v)$Q(v),alt = dist(u) + cost(a);if(alt < dist(v),dist(v) = alt;prev(i,v) = no;prev(u,v) = yes;);););option prev:0:0:7,dist:3:0:6;display"============== Dijkstra (GAMS) ==============",prev,dist;parameter results(*,*);results('Dijkstra (GAMS)','obj') = dist(sink);display results;
- The prev array is here implemented as a 2d set. It is noted that in GAMS, set elements are strings and not numbers.
- Finding a node in \(Q\) with minimum distance from the source is a bit special. We do this in two steps. We are using a singleton set to protect against the case that multiple nodes have the same distance.
- Overall, this implementation did not look too bad.
- The output is listed below.
Listing
---- 97 ============== Dijkstra (GAMS) ============== ---- 97 SET prev i is previous wrt j node1 .node10, node1 .node16, node1 .node24, node1 .node32, node1 .node43, node1 .node45, node2 .node30 node3 .node5 , node3 .node28, node6 .node23, node9 .node26, node9 .node50, node10.node17, node10.node27 node10.node48, node11.node14, node16.node6 , node17.node35, node18.node42, node22.node21, node26.node13 node26.node15, node26.node33, node26.node36, node26.node39, node26.node41, node29.node46, node31.node20 node32.node2 , node32.node4 , node32.node8 , node32.node9 , node32.node11, node32.node18, node32.node22 node32.node29, node32.node38, node32.node47, node38.node3 , node38.node7 , node38.node12, node39.node49 node43.node25, node43.node34, node43.node37, node43.node40, node45.node44, node47.node31, node50.node19 ---- 97 PARAMETER dist distance from source node2 9.465, node3 7.225, node4 10.273, node5 10.690, node6 11.377, node7 9.093 node8 10.352, node9 5.408, node10 4.935, node11 9.075, node12 7.918, node13 9.078 node14 10.576, node15 10.237, node16 4.958, node17 6.860, node18 11.738, node19 9.332 node20 15.809, node21 12.661, node22 9.683, node23 12.467, node24 7.343, node25 8.774 node26 7.534, node27 8.955, node28 9.424, node29 6.406, node30 11.291, node31 10.122 node32 3.319, node33 9.530, node34 10.898, node35 10.377, node36 8.884, node37 11.336 node38 5.416, node39 9.081, node40 14.728, node41 12.327, node42 16.948, node43 6.534 node44 11.847, node45 8.792, node46 9.398, node47 8.475, node48 7.679, node49 10.908 node50 7.803 ---- 101 PARAMETER results obj Dijkstra (GAMS) 10.908
*-----------------------------------------------------* recover path* more complicated than the algorithm itself*-----------------------------------------------------set s /step1*step25/;singleton setsn(i) 'current node'p(i) 'previous node';parameter steps(s,i,j) 'cost for each step taken';option steps:3:0:1;steps(s,i,j) = 0;* initializationn(sink) = yes;* loop from last (sink) to first (source)repeatp(i) = no;* loop finds v given n* body of loop is executed 0 or 1 timesloop(prev(v,n),* make room by moving all steps one downsteps(s,i,j) = steps(s-1,i,j);steps('step1',prev) = cost(prev);p(v) = yes;);n(i) = p(i);until card(n)=0;display steps;
- The set \(s\) is the number of steps of the optimal shortest path. I probably should have made the size of \(s\) larger. A good bound would be the number of nodes.
- The inner loop is not really a loop. It finds, given a single node \(n\), what node \(v\rightarrow n\) is on the optimal shortest path. In general, this is just one node (or no node, if \(n\) is the source node).
- The output is:
---- 134 PARAMETER steps cost for each step taken step1.node1 .node32 3.319 step2.node32.node9 2.089 step3.node9 .node26 2.126 step4.node26.node39 1.547 step5.node39.node49 1.827
Approach 2: Python Implementation
- Python error messages appear out-of-sync in the log file, making it difficult to see what is happening.
- Line numbers in the Python error messages are difficult to interpret. They don't correspond to the line numbers shown in the IDE.
- No Python syntax coloring. All Python code is shown in red (why?).
- No debugging.
- You cannot run the embedded Python code in a standalone Python environment.
- GAMS Python is somewhat barebones, and adding packages is difficult.
- The interfaces between GAMS and Python are not intuitive and are not unified.
- Needless duplications are needed. We have to say 2 times we want to move things back from Python to GAMS: once inside Python and then again in
endEmbeddedCode psteps,pobj
There should be no need to say this twice.
embeddedCode Python: #-- import GAMS data import gams.transfer as gt nodes=list(gams.get('i')) print(f"i={nodes}") source=list(gams.get('source'))[0] print(f"source={source}") sink=list(gams.get('sink'))[0] print(f"sink={sink}") m = gt.Container(gams.db) dfcost = m.data["cost"].records print(dfcost) #-- data structures successors = {v:[] for v in nodes} cost = {} for index, row in dfcost.iterrows(): u = row['i'] v = row['j'] successors[u].append(v) cost[(u,v)] = row['value'] #print(successors) #-- initialization Q=set(nodes) inf=float('inf') dist = {v:inf for v in nodes} dist[source] = 0 prev = {v:None for v in nodes} #-- run algorithm while len(Q)>0: u = min(Q, key=lambda n:dist[n]) Q.remove(u) for v in successors[u]: if v in Q: alt = dist[u] + cost[(u,v)] if alt < dist[v]: dist[v] = alt prev[v] = u print(f"optimal total cost: {dist[sink]}") #-- recover path path = [] u = sink while u is not None: path.insert(0,u) u = prev[u] print(path) #-- move things back to GAMS stepsdata = [(f"step{i}",path[i-1],path[i],cost[(path[i-1],path[i])]) for i in range(1,len(path))] gams.set("psteps",stepsdata) gams.set("pobj",[dist[sink]]) endEmbeddedCode psteps,pobj
- We print the input data to the log so we have a better understanding of how the data arrives. This should look like:
--- Initialize embedded library embpycclib64.dll
--- Execute embedded library embpycclib64.dll
i=['node1', 'node2', 'node3', 'node4', 'node5', 'node6', 'node7', 'node8', 'node9',
source=node1
sink=node49
i j value
0 node1 node10 4.934533
1 node1 node16 4.957717
2 node1 node24 7.343160
3 node1 node32 3.318676
4 node1 node43 6.534425
.. ... ... ...
390 node50 node29 5.084533
391 node50 node34 9.771928
392 node50 node39 8.409906
393 node50 node41 7.046857
394 node50 node47 4.655581
[395 rows x 3 columns] - The set \(i\) arrives as a list of strings. Source and sink are strings. Dfcost is a dataframe (it holds the costs but also the sparsity pattern of the network). We see all 395 arcs arrived safely.
- The first thing to do after importing the data is to create and populate data structures that are better suited for our algorithm. We used dictionaries here so we can keep using strings as nodes. Another approach would have been to map strings to integers and then use integers in the algorithm.
- The optimal solution is:
optimal total cost: 10.907756446999999
['node1', 'node32', 'node9', 'node26', 'node39', 'node49'] - After shipping things back to GAMS, we see in the listing file:
---- 220 ============== Dijkstra (Python) ============== ---- 220 PARAMETER psteps results from python code step1.node1 .node32 3.319 step2.node32.node9 2.089 step3.node9 .node26 2.126 step4.node26.node39 1.547 step5.node39.node49 1.827 ---- 220 PARAMETER pobj = 10.908 objective function value ---- 224 PARAMETER results obj Dijkstra (GAMS) 10.908 Dijkstra (Python) 10.908
Approach 3: Linear Programming Model
| Shortest path LP Model |
|---|
| \[\begin{align}\min& \sum_{(i,j)\in A}\color{darkblue}c_{i,j}\cdot\color{darkred}x_{i,j}\\ & \sum_{j|(j,i)\in A}\color{darkred}x_{j,i} + \color{darkblue}{\mathit{supply}}_i = \sum_{j|(i,j)\in A}\color{darkred}x_{i,j} + \color{darkblue}{\mathit{demand}}_i && \forall i\\ & \color{darkred}x_{i,j} \ge 0 \\ & \color{darkblue}{\mathit{supply}}_i = \begin{cases} 1 & \text{if $i$ is the source node} \\ 0 & \text{otherwise}\end{cases}\\ & \color{darkblue}{\mathit{demand}}_i = \begin{cases} 1 & \text{if $i$ is the sink node} \\ 0 & \text{otherwise}\end{cases} \end{align}\] |
parameters supply(i), demand(i);supply(source) = 1;demand(sink) = 1;positive variables x(i,j) 'flow along i->j';variable z 'objective';equationsobj 'objective'nodebal(i) 'node balance';obj.. z =e= sum((i,j),cost(i,j)*x(i,j));nodebal(i).. sum(a(j,i),x(a)) + supply(i) =e= sum(a(i,j),x(a)) + demand(i);model m /all/;
The results look like:
---- 252 ============== LP model ============== ---- 252 VARIABLE x.L flow along i->j node9 node26 node32 node39 node49 node1 1.000 node9 1.000 node26 1.000 node32 1.000 node39 1.000 ---- 252 VARIABLE z.L = 10.908 objective ---- 255 PARAMETER results obj Dijkstra (GAMS) 10.908 Dijkstra (Python) 10.908 LP 10.908
Again, we want to produce a path from this solution. Not completely trivial, but we can follow the approach from [4].
*-----------------------------------------------------* recover path*-----------------------------------------------------parameterslpsteps(s,i,j) 'results from LP model';singleton set next(i) 'next node';n(i) = source(i);loop(s$card(n),next(i) = x.l(n,i)>0.5;lpsteps(s,n,next) = cost(n,next);n(i) = next(i););option lpsteps:3:0:1;display lpsteps;
---- 273 PARAMETER lpsteps results from LP model step1.node1 .node32 3.319 step2.node32.node9 2.089 step3.node9 .node26 2.126 step4.node26.node39 1.547 step5.node39.node49 1.827
Details
| Dual LP Model |
|---|
| \[\begin{align}\max& \sum_i \left(\color{darkblue}{\mathit{supply}}_i - \color{darkblue}{\mathit{demand}}_i \right)\cdot \color{darkred}\pi_i \\ & \color{darkred}\pi_i -\color{darkred}\pi_j \le \color{darkblue}c_{i,j} && \forall (i,j) \in A \\ & \color{darkred}\pi_{i}\> \textbf{free} \end{align}\] |
Complete GAMS Model
GAMS Model
$onText
Implementation of Dijkstra's algorithm for the shortest path problem. From: https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm
1. Dijkstra's shortest path algorithm implemented in GAMS 2. Dijkstra's shortest path algorithm implemented in Python 3. Shortest path as LP 4. Visualization of network and shortest path $offText
*----------------------------------------------------- * set up network *-----------------------------------------------------
set i 'nodes' /node1*node50/; alias (i,j,v);
set a(i,j) 'arcs of sparse network'; * no self loops and density of 15% a(i,j)$(ord(i)<>ord(j)) = uniform(0,1)<0.15;
singleton sets source(i) /node1/ sink(i) /node49/ ; * sink is chosen so optimal path needs quite a few hops
parameter cost(i,j) 'random costs (lengths)'; cost(a) = uniform(1,10);
option a:0:0:7, cost:3:0:6; display$(card(i)<=50) "============== data ==============", i,source,sink,a,cost;
*===================================================== * 1. Dijkstra GAMS algorithm *=====================================================
*----------------------------------------------------- * data structures *-----------------------------------------------------
sets prev(i,j) 'i is previous wrt j' Q(i) 'unvisited nodes' ; parameter dist(i) 'distance from source';
* initialize Q(i) = yes; dist(i) = INF; dist(source) = 0;
*----------------------------------------------------- * algorithm *-----------------------------------------------------
scalar mindist; singleton set u(i); scalar alt;
* if multiple assignments to singleton set, * only the last one will remain. This will * guarantee a singleton set will hold * 0 or 1 elements. option strictSingleton=0;
while(card(Q)>0,
* find node in Q with minimum distance from the source mindist = smin(q,dist(q)); u(Q) = dist(Q)=mindist;
* remove u from Q Q(u) = no;
* loop over neighbors of u in Q loop(a(u,v)$Q(v), alt = dist(u) + cost(a); if(alt < dist(v), dist(v) = alt; prev(i,v) = no; prev(u,v) = yes; ); ); );
option prev:0:0:7,dist:3:0:6; display "============== Dijkstra (GAMS) ==============", prev,dist; parameter results(*,*); results('Dijkstra (GAMS)','obj') = dist(sink); display results;
*----------------------------------------------------- * recover path * more complicated than the algorithm itself *-----------------------------------------------------
set s /step1*step25/; singleton sets n(i) 'current node' p(i) 'previous node' ; parameter steps(s,i,j) 'cost for each step taken'; option steps:3:0:1; steps(s,i,j) = 0;
* initialization n(sink) = yes;
* loop from last (sink) to first (source) repeat p(i) = no; * loop finds v given n * body of loop is executed 0 or 1 times loop(prev(v,n), * make room by moving all steps one down steps(s,i,j) = steps(s-1,i,j); steps('step1',prev) = cost(prev); p(v) = yes; ); n(i) = p(i); until card(n)=0;
display steps;
*===================================================== * 2. Dijkstra Python algorithm *=====================================================
parameters psteps(s,i,j) 'results from python code' pobj 'objective function value' ;
* just in case we have arcs with costs=0 cost(a)$(cost(a)=0) = eps;
embeddedCode Python:
#-- import GAMS data import gams.transfer as gt
nodes=list(gams.get('i')) print(f"i={nodes}")
source=list(gams.get('source'))[0] print(f"source={source}")
sink=list(gams.get('sink'))[0] print(f"sink={sink}")
m = gt.Container(gams.db) dfcost = m.data["cost"].records print(dfcost)
#-- data structures successors = {v:[] for v in nodes} cost = {} for index, row in dfcost.iterrows(): u = row['i'] v = row['j'] successors[u].append(v) cost[(u,v)] = row['value']
#print(successors)
#-- initialization Q=set(nodes) inf=float('inf') dist = {v:inf for v in nodes} dist[source] = 0 prev = {v:None for v in nodes}
#-- run algorithm while len(Q)>0: u = min(Q, key=lambda n:dist[n]) Q.remove(u) for v in successors[u]: if v in Q: alt = dist[u] + cost[(u,v)] if alt < dist[v]: dist[v] = alt prev[v] = u print(f"optimal total cost: {dist[sink]}")
#-- recover path path = [] u = sink while u is not None: path.insert(0,u) u = prev[u] print(path)
#-- move things back to GAMS
stepsdata = [(f"step{i}",path[i-1],path[i],cost[(path[i-1],path[i])]) for i in range(1,len(path))] gams.set("psteps",stepsdata) gams.set("pobj",[dist[sink]])
endEmbeddedCode psteps,pobj
option psteps:3:0:1; display "============== Dijkstra (Python) ==============", psteps,pobj;
results('Dijkstra (Python)','obj') = pobj; display results;
*===================================================== * 3. LP Model *=====================================================
parameters supply(i), demand(i); supply(source) = 1; demand(sink) = 1;
positive variables x(i,j) 'flow along i->j'; variable z 'objective';
equations obj 'objective' nodebal(i) 'node balance' ;
obj.. z =e= sum((i,j),cost(i,j)*x(i,j)); nodebal(i).. sum(a(j,i),x(a)) + supply(i) =e= sum(a(i,j),x(a)) + demand(i);
model m /all/; solve m minimizing z using lp;
display "============== LP model ==============", x.l,z.l; results('LP','obj') = z.l; display results;
*----------------------------------------------------- * recover path *-----------------------------------------------------
parameters lpsteps(s,i,j) 'results from LP model' ; singleton set next(i) 'next node';
n(i) = source(i); loop(s$card(n), next(i) = x.l(n,i)>0.5; lpsteps(s,n,next) = cost(n,next); n(i) = next(i); ); option lpsteps:3:0:1; display lpsteps;
*===================================================== * 4. Visualization *=====================================================
*----------------------------------------------------------- * Plot network using put statement * We generate some JavaScript code to hold the data * For documentation on plotting library, see: * https://js.cytoscape.org/ *-----------------------------------------------------------
$set htmlfile network.html $set datafile networkdata.js
file fln /%datafile%/; put fln; put 'networkdata=['/; loop(i, put$(ord(i)>1) ","; put "{data:{id:",(ord(i)):0:0; put$(supply(i)=0 and demand(i)=0) ",color:'black'"; put$(supply(i)>0) ",color:'blue'"; put$(demand(i)>0) ",color:'green'"; put$(supply(i)=0 and demand(i)=0) ",size:1"; put$(supply(i)>0 or demand(i)>0) ",size:2"; put "}}"/; ); loop(a(i,j), put ",{data:{id:'",(ord(i)):0:0,'-',(ord(j)):0:0,"',"; put "source:",(ord(i)):0:0,",target:",(ord(j)):0:0; put$(x.l(i,j)>0.5) ",color:'red',width:0.2"; put$(x.l(i,j)<0.5) ",color:'grey',width:0.1"; put "}}"/; ); put '];'/;
put "table='<h4>GAMS Shortest Path</h4>" "Number of nodes: ",(card(i)):0:0,"<br>" "Number of arcs: ",(card(a)):0:0,"<br><br>" "<table><tr><th>From</th><th>To</th><th>cost</th></tr>'+"/; loop((s,i,j)$lpsteps(s,i,j), put "'<tr><td>",i.tl:0,"</td><td>",j.tl:0,"</td><td><pre>",cost(i,j):10:3,"</pre></td></tr>'+"/; ); put "'<tr><td colspan=2>Total cost</td><td><pre>",z.l:10:3,"</pre></td></tr>'+"/;
put "'</table><br>';"; putclose;
$ontext
With an internet connection you can use: <script src="https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.30.4/cytoscape.min.js"></script> Without use a local .js file: <script src="cytoscape.min.js"></script>
$offtext
$onecho > %htmlfile% <html> <head> <title>GAMS Shortest Path Network</title> <script src="https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.30.4/cytoscape.min.js"></script> <script src="networkdata.js"></script> </head>
<style> #cy { width: calc(100% - 200px); height: 100%; position: absolute; top: 0px; left: 200px; }
table,th, td { border: 1px solid black; border-collapse: collapse; padding-left: 10px; padding-right: 10px; }
</style>
<body> <div id="cy"></div> <div id="my-table"></div> <div><button onclick="buttonClicked()" id="btn">Animate Flow</button></div> <script> document.getElementById('my-table').innerHTML = table var cy = cytoscape({ container: document.getElementById('cy'), elements: networkdata, style: [ { selector: 'node', style: { 'background-color': 'data(color)', label: 'data(id)', width: 'data(size)', height: 'data(size)', 'font-size': 1.5 } }, { selector: 'edge', style: { 'width': 'data(width)', 'line-color': 'data(color)', 'mid-target-arrow-shape': 'triangle', 'mid-target-arrow-color': 'data(color)', 'arrow-scale': 0.15, 'curve-style': 'bezier' }
} ], layout: { name: 'cose' } });
const loopAnimation = (ele, i) => { const offset = { style: {'line-dash-offset': -100 * i } }; const duration = { duration: 10000 }; return ele.animation(offset, duration).play() .promise('complete') .then(() => loopAnimation(ele, i + 1)); };
var reds = cy.edges().filter(function(ele) { return ele.data('color') == 'red'; }); reds.forEach((edge) => { loopAnimation(edge, 1); });
var animated = false; const btn = document.getElementById("btn");
function buttonClicked(b) { if (animated) { reds.style({'line-style':'solid', 'width':0.2}); btn.innerHTML = "Animate flow"; } else { reds.style({'line-style':'dashed', 'line-dash-pattern': [0.6, 1.1], 'width':0.4}); btn.innerHTML = "Stop animation"; } animated = !animated; } </script> </body> </html> $offecho
execute "%htmlfile%"; |
References
- E. W. Dijkstra, A note on two problems in connexion with graphs., Numer. Math. 1, 269–271 (1959).
- Edsger W. Dijkstra, https://en.wikipedia.org/wiki/Edsger_W._Dijkstra
- Dijkstra's Algorithm, https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm
- A network model: Pyomo vs GAMS, https://yetanothermathprogrammingconsultant.blogspot.com/2021/08/a-network-model-pyomo-vs-gams.html
- https://www.modelworks.ch/dijkstra.htm, another implementation of Dijkstra's algorithm in GAMS.
 |
| Click to enlarge |
No comments:
Post a Comment