Monday, February 24, 2025

Programming vs Modeling

I found another interesting GAMS model. Here, I wanted to demonstrate the difference between programming (implementing algorithms in a programming language) and modeling (building models and implementing them in a modeling tool). The problem is solving a shortest path problem given a sparse network with random costs. We compare three strategies:

  1. Implement Dijkstra's shortest path algorithm in GAMS,
  2. implement Dijkstra's shortest path algorithm in Python and
  3. use a simple LP model.  

Data


For this exercise, I generate a random, sparse directed graph, with some random costs (lengths) for the arcs. The GAMS code looks like:

*-----------------------------------------------------
* 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;


A directed graph is mathematically described as \(G=(V,E)\). We have a set of vertices (node) \(V\) and a set of directed edges (arcs) \(E\). Here, I used \(i\) to denote our set of nodes and \(a\) to describe our arcs. We have 50 nodes, and we generated 395 arcs. We don't allow self-loops (that is an arc \(i \rightarrow i\)). We aimed for 15% density, and we ended up with about 16%. 

Notes:
  • Singleton sets can be used to enforce that they only contain 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. 

The data set looks a bit more impressive when we print it out.

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



Dijkstra's algorithm was first published in a very short paper [1] in 1959.  


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".
I used the pseudo code from [3]:


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


GAMS is a modeling language. We can, however, write little algorithms in GAMS. Obviously, one would not write full-blown applications (like a word processor or an LP solver) in a language like GAMS. But Dijkstra's algorithm may be quite doable. Here is my attempt.

*-----------------------------------------------------
* 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;

Notes:
  • 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


The code to recover the path from source to sink is a bit more complicated. Actually, more complicated than the Dijkstra algorithm itself.

*-----------------------------------------------------
* 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;

Notes:
  • 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

The optimal objective value is the sum of the costs (lengths) of each step. The extra set \(s\) is needed here to maintain the correct ordering of the \((i,j)\) pairs.

Conclusion: GAMS is not a real programming language, but it is expressive enough to be able to implement a (simple version) of Dijkstra's algorithm.

Approach 2: Python Implementation


GAMS allows you to write pieces of Python code inside a GAMS model. This is the so-called Embedded Python facility. In practice, it is not that easy to use. Some reasons:
  • 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.
Because of all these drawbacks, it is just not very pleasant (or fun) to use. 

Here is the code (with added syntax coloring):

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

Again, we follow the algorithm from [3].

Notes:
  • 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

Reassuringly, we arrived at the same optimal solution.

Conclusion: the Python code is a bit more compact than the GAMS implementation. However, there is quite some extra work involved in importing, exporting, and transforming the GAMS data.

Approach 3: Linear Programming Model


In the previous two sections, we implemented Dijkstra's algorithm in GAMS and Python. Here, we use as an alternative an LP formulation. This is very different. The model looks like:

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}\]

The corresponding GAMS code is straightforward once we have this mathematical 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/;

Instead of writing down how to solve the problem, here we write down what to solve. It is a very different way of looking at the problem.

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
*-----------------------------------------------------

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;


This gives us:

----    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

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
  1. E. W. Dijkstra,  A note on two problems in connexion with graphs., Numer. Math. 1, 269–271 (1959).
  2. Edsger W. Dijkstra, https://en.wikipedia.org/wiki/Edsger_W._Dijkstra
  3. Dijkstra's Algorithm, https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm
  4. A network model: Pyomo vs GAMS, https://yetanothermathprogrammingconsultant.blogspot.com/2021/08/a-network-model-pyomo-vs-gams.html
  5. https://www.modelworks.ch/dijkstra.htm, another implementation of Dijkstra's algorithm in GAMS.

Click to enlarge

Note that the locations of the nodes are arbitrary in this picture and the lengths of the arcs are not related to the cost vector that is used in the optimization model.
Posted by at
Labels: , ,

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)