1# Copyright 2011 Hakan Kjellerstrand hakank@gmail.com
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7#     http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
14"""
15
16  3 jugs problem using MIP in Google or-tools.
17
18  A.k.a. water jugs problem.
19
20  Problem from Taha 'Introduction to Operations Research',
21  page 245f .
22
23  Compare with the CP model:
24     http://www.hakank.org/google_or_tools/3_jugs_regular
25
26  This model was created by Hakan Kjellerstrand (hakank@gmail.com)
27  Also see my other Google CP Solver models:
28  http://www.hakank.org/google_or_tools/
29"""
30import sys
31from ortools.linear_solver import pywraplp
32
33
34def main(sol='CBC'):
35
36  # Create the solver.
37
38  print('Solver: ', sol)
39
40  # using GLPK
41  if sol == 'GLPK':
42    solver = pywraplp.Solver('CoinsGridGLPK',
43                             pywraplp.Solver.GLPK_MIXED_INTEGER_PROGRAMMING)
44  else:
45    # Using CBC
46    solver = pywraplp.Solver('CoinsGridCBC',
47                             pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)
48
49  #
50  # data
51  #
52  n = 15
53  start = 0  # start node
54  end = 14  # end node
55  M = 999  # a large number
56
57  nodes = [
58      '8,0,0',  # start
59      '5,0,3',
60      '5,3,0',
61      '2,3,3',
62      '2,5,1',
63      '7,0,1',
64      '7,1,0',
65      '4,1,3',
66      '3,5,0',
67      '3,2,3',
68      '6,2,0',
69      '6,0,2',
70      '1,5,2',
71      '1,4,3',
72      '4,4,0'  # goal!
73  ]
74
75  # distance
76  d = [[M, 1, M, M, M, M, M, M, 1, M, M, M, M, M, M],
77       [M, M, 1, M, M, M, M, M, M, M, M, M, M, M, M],
78       [M, M, M, 1, M, M, M, M, 1, M, M, M, M, M, M],
79       [M, M, M, M, 1, M, M, M, M, M, M, M, M, M, M],
80       [M, M, M, M, M, 1, M, M, 1, M, M, M, M, M, M],
81       [M, M, M, M, M, M, 1, M, M, M, M, M, M, M, M],
82       [M, M, M, M, M, M, M, 1, 1, M, M, M, M, M, M],
83       [M, M, M, M, M, M, M, M, M, M, M, M, M, M, 1],
84       [M, M, M, M, M, M, M, M, M, 1, M, M, M, M, M],
85       [M, 1, M, M, M, M, M, M, M, M, 1, M, M, M, M],
86       [M, M, M, M, M, M, M, M, M, M, M, 1, M, M, M],
87       [M, 1, M, M, M, M, M, M, M, M, M, M, 1, M, M],
88       [M, M, M, M, M, M, M, M, M, M, M, M, M, 1, M],
89       [M, 1, M, M, M, M, M, M, M, M, M, M, M, M, 1],
90       [M, M, M, M, M, M, M, M, M, M, M, M, M, M, M]]
91
92  #
93  # variables
94  #
95
96  # requirements (right hand statement)
97  rhs = [solver.IntVar(-1, 1, 'rhs[%i]' % i) for i in range(n)]
98
99  x = {}
100  for i in range(n):
101    for j in range(n):
102      x[i, j] = solver.IntVar(0, 1, 'x[%i,%i]' % (i, j))
103
104  out_flow = [solver.IntVar(0, 1, 'out_flow[%i]' % i) for i in range(n)]
105  in_flow = [solver.IntVar(0, 1, 'in_flow[%i]' % i) for i in range(n)]
106
107  # length of path, to be minimized
108  z = solver.Sum(
109      [d[i][j] * x[i, j] for i in range(n) for j in range(n) if d[i][j] < M])
110
111  #
112  # constraints
113  #
114
115  for i in range(n):
116    if i == start:
117      solver.Add(rhs[i] == 1)
118    elif i == end:
119      solver.Add(rhs[i] == -1)
120    else:
121      solver.Add(rhs[i] == 0)
122
123  # outflow constraint
124  for i in range(n):
125    solver.Add(
126        out_flow[i] == solver.Sum([x[i, j] for j in range(n) if d[i][j] < M]))
127
128  # inflow constraint
129  for j in range(n):
130    solver.Add(
131        in_flow[j] == solver.Sum([x[i, j] for i in range(n) if d[i][j] < M]))
132
133  # inflow = outflow
134  for i in range(n):
135    solver.Add(out_flow[i] - in_flow[i] == rhs[i])
136
137  # objective
138  objective = solver.Minimize(z)
139
140  #
141  # solution and search
142  #
143  solver.Solve()
144
145  print()
146  print('z: ', int(solver.Objective().Value()))
147
148  t = start
149  while t != end:
150    print(nodes[t], '->', end=' ')
151    for j in range(n):
152      if x[t, j].SolutionValue() == 1:
153        print(nodes[j])
154        t = j
155        break
156
157  print()
158  print('walltime  :', solver.WallTime(), 'ms')
159  if sol == 'CBC':
160    print('iterations:', solver.Iterations())
161
162
163if __name__ == '__main__':
164
165  sol = 'CBC'
166  if len(sys.argv) > 1:
167    sol = sys.argv[1]
168    if sol != 'GLPK' and sol != 'CBC':
169      print('Solver must be either GLPK or CBC')
170      sys.exit(1)
171
172  main(sol)
173