Solving a large system of equations for the smallest integer values

I’m trying to find the smallest possible non-negative integer solutions to hundreds of fairly large equations systems . An average example would be the following:

expr = {x51 + x61 + y41 + y51 + y61 >= 1,
 x51 + x61 + y41 + y51 + y61 == x12 + x13 + x14 + x15 + x16 + y13 + y14 + y15 + y16 == x12 + x23 + x34 + x45 + x56 + y13 + y24 + y35 + y46,
 x12 + x62 + y52 + y62 >= 1,
 x12 + x62 + y52 + y62 == x23 + x24 + x25 + x26 + y24 + y25 + y26 == x13 + x24 + x35 + x46 + y14 + y25 + y36,
 x13 + x23 + y13 + y63 >= 1,
 x13 + x23 + y13 + y63 == x34 + x35 + x36 + y35 + y36 == x14 + x25 + x36 + y15 + y26,
 x14 + x24 + x34 + y14 + y24 >= 1,
 x14 + x24 + x34 + y14 + y24 == x45 + x46 + y41 + y46 == x15 + x26 + y16 + y61,
 x15 + x25 + x35 + x45 + y15 + y25 + y35 >= 1,
 x15 + x25 + x35 + x45 + y15 + y25 + y35 == x51 + x56 + y51 + y52 == x16 + x61 + y51 + y62,
 x16 + x26 + x36 + x46 + x56 + y16 + y26 + y36 + y46 >= 1,
 x16 + x26 + x36 + x46 + x56 + y16 + y26 + y36 + y46 == x61 + x62 + y61 + y62 + y63 == x51 + x62 + y41 + y52 + y63,
 x51 + x61 + x62 == y13 + y14 + y15 + y16 + y24 + y25 + y26 + y35 + y36 + y46,
 y41 + y51 + y52 + y61 + y62 + y63 >= 1 [Implies] x51 + x61 + x62 >= 1,
 {x12, x13, x14, x15, x16, x23, x24, x25, x26, x34, x35, x36, x45, x46, x51, x56, x61, x62, y13, y14, y15, y16, y24, y25, y26, y35, y36, y41, y46, y51, y52, y61, y62, y63} [Element] NonNegativeIntegers}; 

Running Solve[] for all the variables does not terminate after several hours. On smaller systems than the one above, I was able to use Normal[Solve[expr,vars]] which gave me a set of rules as the one below:

sols = {{a -> 1 + 3 C[1] + 2 C[2] + 2 C[3] + C[4] + 2 C[5] + C[6] + 3 C[7] + C[8] + C[9], 
 b -> C[1] + 3 C[2] + 3 C[3] + C[4] + 3 C[5] + C[6], 
 c -> 2 C[1] + C[2] + 2 C[3] + C[4] + 2 C[7] + C[8], 
 d -> 1 + 2 C[1] + C[2] + 2 C[5] + C[6] + 2 C[7] + C[9], 
 e -> 3 + C[4] + C[6] + 2 C[7] + 3 C[8] + 3 C[9], 
 m -> C[1] + C[7], 
 n -> C[1] + C[7], 
 o -> 2 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 p -> 1 + C[2] + 2 C[5] + C[6] + C[9], 
 q -> 2 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 r -> C[2] + 2 C[3] + C[4] + C[8]}, 
 {a -> 1 + 3 C[1] + 2 C[2] + 2 C[3] + C[4] + 2 C[5] + C[6] + 3 C[7] + C[ 8] + C[9], 
 b -> C[1] + 3 C[2] + 3 C[3] + C[4] + 3 C[5] + C[6], 
 c -> 1 + 2 C[1] + C[2] + 2 C[3] + C[4] + 2 C[7] + C[8], 
 d -> 2 C[1] + C[2] + 2 C[5] + C[6] + 2 C[7] + C[9], 
 e -> 3 + C[4] + C[6] + 2 C[7] + 3 C[8] + 3 C[9], 
 m -> C[1] + C[7], 
 n -> C[1] + C[7], 
 o -> 2 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 p -> C[2] + 2 C[5] + C[6] + C[9], 
 q -> 2 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 r -> 1 + C[2] + 2 C[3] + C[4] + C[8]}, 
 {a -> 1 + 3 C[1] + 2 C[2] + 2 C[3] + C[4] + 2 C[5] + C[6] + 3 C[7] + C[ 8] + C[9], 
 b -> 1 + C[1] + 3 C[2] + 3 C[3] + C[4] + 3 C[5] + C[6], 
 c -> 2 C[1] + C[2] + 2 C[3] + C[4] + 2 C[7] + C[8], 
 d -> 1 + 2 C[1] + C[2] + 2 C[5] + C[6] + 2 C[7] + C[9], 
 e -> 1 + C[4] + C[6] + 2 C[7] + 3 C[8] + 3 C[9], 
 m -> C[1] + C[7], 
 n -> C[1] + C[7], 
 o -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 p -> 1 + C[2] + 2 C[5] + C[6] + C[9], 
 q -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 r -> C[2] + 2 C[3] + C[4] + C[8]}, 
 {a -> 1 + 3 C[1] + 2 C[2] + 2 C[3] + C[4] + 2 C[5] + C[6] + 3 C[7] + C[ 8] + C[9], 
 b -> 1 + C[1] + 3 C[2] + 3 C[3] + C[4] + 3 C[5] + C[6], 
 c -> 1 + 2 C[1] + C[2] + 2 C[3] + C[4] + 2 C[7] + C[8], 
 d -> 2 C[1] + C[2] + 2 C[5] + C[6] + 2 C[7] + C[9], 
 e -> 1 + C[4] + C[6] + 2 C[7] + 3 C[8] + 3 C[9], 
 m -> C[1] + C[7], 
 n -> C[1] + C[7], 
 o -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 p -> C[2] + 2 C[5] + C[6] + C[9], 
 q -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 r -> 1 + C[2] + 2 C[3] + C[4] + C[8]}, 
 {a -> 2 + 3 C[1] + 2 C[2] + 2 C[3] + C[4] + 2 C[5] + C[6] + 3 C[7] + C[ 8] + C[9], 
 b -> 3 + C[1] + 3 C[2] + 3 C[3] + C[4] + 3 C[5] + C[6], 
 c -> 2 C[1] + C[2] + 2 C[3] + C[4] + 2 C[7] + C[8], 
 d -> 2 + 2 C[1] + C[2] + 2 C[5] + C[6] + 2 C[7] + C[9], 
 e -> C[4] + C[6] + 2 C[7] + 3 C[8] + 3 C[9], 
 m -> C[1] + C[7], 
 n -> C[1] + C[7], 
 o -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 p -> 2 + C[2] + 2 C[5] + C[6] + C[9], 
 q -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 r -> C[2] + 2 C[3] + C[4] + C[8]}, 
 {a -> 2 + 3 C[1] + 2 C[2] + 2 C[3] + C[4] + 2 C[5] + C[6] + 3 C[7] + C[ 8] + C[9], 
 b -> 3 + C[1] + 3 C[2] + 3 C[3] + C[4] + 3 C[5] + C[6], 
 c -> 1 + 2 C[1] + C[2] + 2 C[3] + C[4] + 2 C[7] + C[8], 
 d -> 1 + 2 C[1] + C[2] + 2 C[5] + C[6] + 2 C[7] + C[9], 
 e -> C[4] + C[6] + 2 C[7] + 3 C[8] + 3 C[9], 
 m -> C[1] + C[7], 
 n -> C[1] + C[7], 
 o -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 p -> 1 + C[2] + 2 C[5] + C[6] + C[9], 
 q -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 r -> 1 + C[2] + 2 C[3] + C[4] + C[8]}, 
 {a -> 2 + 3 C[1] + 2 C[2] + 2 C[3] + C[4] + 2 C[5] + C[6] + 3 C[7] + C[ 8] + C[9], 
 b -> 3 + C[1] + 3 C[2] + 3 C[3] + C[4] + 3 C[5] + C[6], 
 c -> 2 + 2 C[1] + C[2] + 2 C[3] + C[4] + 2 C[7] + C[8], 
 d -> 2 C[1] + C[2] + 2 C[5] + C[6] + 2 C[7] + C[9], 
 e -> C[4] + C[6] + 2 C[7] + 3 C[8] + 3 C[9], 
 m -> C[1] + C[7], 
 n -> C[1] + C[7], 
 o -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 p -> C[2] + 2 C[5] + C[6] + C[9], 
 q -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 r -> 2 + C[2] + 2 C[3] + C[4] + C[8]}, 
 {a -> 3 + 3 C[1] + 2 C[2] + 2 C[3] + C[4] + 2 C[5] + C[6] + 3 C[7] + C[ 8] + C[9], 
 b -> C[1] + 3 C[2] + 3 C[3] + C[4] + 3 C[5] + C[6], 
 c -> 2 + 2 C[1] + C[2] + 2 C[3] + C[4] + 2 C[7] + C[8], 
 d -> 2 + 2 C[1] + C[2] + 2 C[5] + C[6] + 2 C[7] + C[9], 
 e -> 2 + C[4] + C[6] + 2 C[7] + 3 C[8] + 3 C[9], 
 m -> 1 + C[1] + C[7], 
 n -> 1 + C[1] + C[7], 
 o -> 2 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 p -> C[2] + 2 C[5] + C[6] + C[9], 
 q -> 2 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 r -> C[2] + 2 C[3] + C[4] + C[8]}, 
 {a -> 3 + 3 C[1] + 2 C[2] + 2 C[3] + C[4] + 2 C[5] + C[6] + 3 C[7] + C[ 8] + C[9], 
 b -> 1 + C[1] + 3 C[2] + 3 C[3] + C[4] + 3 C[5] + C[6], 
 c -> 2 + 2 C[1] + C[2] + 2 C[3] + C[4] + 2 C[7] + C[8], 
 d -> 2 + 2 C[1] + C[2] + 2 C[5] + C[6] + 2 C[7] + C[9], 
 e -> C[4] + C[6] + 2 C[7] + 3 C[8] + 3 C[9], 
 m -> 1 + C[1] + C[7], 
 n -> 1 + C[1] + C[7], 
 o -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 p -> C[2] + 2 C[5] + C[6] + C[9], 
 q -> 1 + C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + 2 C[7] + 2 C[8] + 2 C[9], 
 r -> C[2] + 2 C[3] + C[4] + C[8]}} 

I guessed the smallest solutions would be those where all the C[x]’s would be equal to zero, so I created the following array…

czeros = Reap[For[i = 1, i < 60, i++, Sow[C[i] -> 0]]]

… and substituted it back in the set of rules that resulted from the Normal[Solve[]]. I then selected the result that gave the smallest sum of all the variables:

varsum = a + b + c + d + e + m + n + o + p + q + r
list = varsum /. sols /. czeros
min = If[Head[Position[list, Min[list]]] == List, Position[list, Min[list]][[1, 1]], Position[list, Min[list]]]
sols[[min]] /. czeros

Having set the scenario, I have two question:

(1) Does the method above indeed gives the solution with the smallest non-negative sum of the variables?

(2) Is there any alternative method that could be used on systems as large as the first example and that would give results in a reasonable amount of time?

It’s my first time using Mathematica, so I’m aware that the code that I’m using might be horribly ineficcient. I have done my best trying to search the documentation for alternatives, but haven´t been able to find any. Thanks in advance for any help.

constraints = {x51 + x61 + y41 + y51 + y61 >= 1, 
   x51 + x61 + y41 + y51 + y61 == 
    x12 + x13 + x14 + x15 + x16 + y13 + y14 + y15 + y16 == 
    x12 + x23 + x34 + x45 + x56 + y13 + y24 + y35 + y46, 
   x12 + x62 + y52 + y62 >= 1, 
   x12 + x62 + y52 + y62 == x23 + x24 + x25 + x26 + y24 + y25 + y26 ==
     x13 + x24 + x35 + x46 + y14 + y25 + y36, 
   x13 + x23 + y13 + y63 >= 1, 
   x13 + x23 + y13 + y63 == x34 + x35 + x36 + y35 + y36 == 
    x14 + x25 + x36 + y15 + y26, x14 + x24 + x34 + y14 + y24 >= 1, 
   x14 + x24 + x34 + y14 + y24 == x45 + x46 + y41 + y46 == 
    x15 + x26 + y16 + y61, 
   x15 + x25 + x35 + x45 + y15 + y25 + y35 >= 1, 
   x15 + x25 + x35 + x45 + y15 + y25 + y35 == x51 + x56 + y51 + y52 ==
     x16 + x61 + y51 + y62, 
   x16 + x26 + x36 + x46 + x56 + y16 + y26 + y36 + y46 >= 1, 
   x16 + x26 + x36 + x46 + x56 + y16 + y26 + y36 + y46 == 
    x61 + x62 + y61 + y62 + y63 == x51 + x62 + y41 + y52 + y63, 
   x51 + x61 + x62 == 
    y13 + y14 + y15 + y16 + y24 + y25 + y26 + y35 + y36 + y46, 
   y41 + y51 + y52 + y61 + y62 + y63 >= 1 [Implies] 
    x51 + x61 + x62 >= 1
   };

variables = {x12, x13, x14, x15, x16, x23, x24, x25, x26, x34, x35, 
   x36, x45, x46, x51, x56, x61, x62, y13, y14, y15, y16, y24, y25, 
   y26, y35, y36, y41, y46, y51, y52, y61, y62, y63};

{minvarsum, solution} = Minimize[{Total[variables], And @@ constraints}, 
  variables, NonNegativeIntegers]

(** result: 
  {7, {x12 -> 0, x13 -> 0, x14 -> 0, x15 -> 0, x16 -> 0, x23 -> 0, 
  x24 -> 0, x25 -> 0, x26 -> 1, x34 -> 0, x35 -> 0, x36 -> 1, 
  x45 -> 1, x46 -> 0, x51 -> 1, x56 -> 0, x61 -> 0, x62 -> 0, 
  y13 -> 0, y14 -> 1, y15 -> 0, y16 -> 0, y24 -> 0, y25 -> 0, 
  y26 -> 0, y35 -> 0, y36 -> 0, y41 -> 0, y46 -> 0, y51 -> 0, 
  y52 -> 0, y61 -> 0, y62 -> 1, y63 -> 1}}
**)