How to solve flow game using Google OR tools?

As I don't have experience with OR-tools, here is an approach with Z3.

• The initial board is represented by numbers for the end points, one number for each color. The idea is a bit similar to how Sudoku is represented.
• Each other cell on the board will get either a value for zero, or a number. This number should be equal to exactly two of its neighbors.
• The initial endpoints should have exactly one neighbor with its color.

from z3 import Solver, Sum, Int, If, And, Or, sat

def plot_solution(S):
    import matplotlib.pyplot as plt

    ax = plt.gca()
    colors = plt.cm.tab10.colors
    for i in range(M):
        for j in range(N):
            if board[i][j] != 0:
                ax.scatter(j, i, s=500, color=colors[board[i][j]])
            if S[i][j] != 0:
                for k in range(M):
                    for l in range(N):
                        if abs(k - i) + abs(l - j) == 1 and S[i][j] == S[k][l]:
                            ax.plot([j, l], [i, k], color=colors[S[i][j]], lw=15)
    ax.set_ylim(M - 0.5, -0.5)
    ax.set_xlim(-0.5, N - 0.5)
    ax.set_aspect('equal')
    ax.set_facecolor('black')
    ax.set_yticks([i + 0.5 for i in range(M - 1)], minor=True)
    ax.set_xticks([j + 0.5 for j in range(N - 1)], minor=True)
    ax.grid(b=True, which='minor', color='white')
    ax.set_xticks([])
    ax.set_yticks([])
    ax.tick_params(axis='both', which='both', length=0)
    plt.show()

board = [[1, 0, 0, 2, 3],
         [0, 0, 0, 4, 0],
         [0, 0, 4, 0, 0],
         [0, 2, 3, 0, 5],
         [0, 1, 5, 0, 0]]
M = len(board)
N = len(board[0])
B = [[Int(f'B_{i}_{j}') for j in range(N)] for i in range(M)]
s = Solver()
s.add(([If(board[i][j] != 0, B[i][j] == board[i][j], And(B[i][j] >= 0, B[i][j] < 10))
        for j in range(N) for i in range(M)]))
for i in range(M):
    for j in range(N):
        same_neighs_ij = Sum([If(B[i][j] == B[k][l], 1, 0)
                              for k in range(M) for l in range(N) if abs(k - i) + abs(l - j) == 1])
        if board[i][j] != 0:
            s.add(same_neighs_ij == 1)
        else:
            s.add(Or(same_neighs_ij == 2, B[i][j] == 0))

if s.check() == sat:
    m = s.model()
    S = [[m[B[i][j]].as_long() for j in range(N)] for i in range(M)]
    print(S)
    plot_solution(S)

Solution:

[[1, 2, 2, 2, 3],
 [1, 2, 4, 4, 3],
 [1, 2, 4, 3, 3],
 [1, 2, 3, 3, 5],
 [1, 1, 5, 5, 5]]

As mentioned in the comments, a possible requirement is that all cells would need to be colored. This would need a more complicated approach. Here is an example of such a configuration for which the above code could create a solution that connects all end points without touching all cells:

board = [[0, 1, 2, 0, 0, 0, 0],
         [1, 3, 4, 0, 3, 5, 0],
         [0, 0, 0, 0, 0, 0, 0],
         [0, 2, 0, 4, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 0],
         [0, 5, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 0]]

The best way is probably with AddCircuit. This constraint takes a directed graph where each arc is labeled with a literal and requires that the arcs labeled true form a subgraph where each node has in- and out-degree 1, and further that there is at most one cycle that is not a self-loop. By forcing an arc from the end to the beginning, we can use this constraint type to require that there is a single path from the beginning to the end.

The documentation is somewhat poor, so here's a working code sample. I'll leave the drawing part to you.

import collections
from ortools.sat.python import cp_model


def validate_board_and_count_colors(board):
    assert isinstance(board, list)
    assert all(isinstance(row, list) for row in board)
    assert len(set(map(len, board))) == 1
    colors = collections.Counter(square for row in board for square in row)
    del colors[0]
    assert all(count == 2 for count in colors.values())
    num_colors = len(colors)
    assert set(colors.keys()) == set(range(1, num_colors + 1))
    return num_colors


def main(board):
    num_colors = validate_board_and_count_colors(board)
    model = cp_model.CpModel()
    solution = [
        [square or model.NewIntVar(1, num_colors, "") for (j, square) in enumerate(row)]
        for (i, row) in enumerate(board)
    ]
    true = model.NewBoolVar("")
    model.AddBoolOr([true])
    for color in range(1, num_colors + 1):
        endpoints = []
        arcs = []
        for i, row in enumerate(board):
            for j, square in enumerate(row):
                if square == color:
                    endpoints.append((i, j))
                else:
                    arcs.append(((i, j), (i, j)))
                if i < len(board) - 1:
                    arcs.append(((i, j), (i + 1, j)))
                if j < len(row) - 1:
                    arcs.append(((i, j), (i, j + 1)))
        (i1, j1), (i2, j2) = endpoints
        k1 = i1 * len(row) + j1
        k2 = i2 * len(row) + j2
        arc_variables = [(k2, k1, true)]
        for (i1, j1), (i2, j2) in arcs:
            k1 = i1 * len(row) + j1
            k2 = i2 * len(row) + j2
            edge = model.NewBoolVar("")
            if k1 == k2:
                model.Add(solution[i1][j1] != color).OnlyEnforceIf(edge)
                arc_variables.append((k1, k1, edge))
            else:
                model.Add(solution[i1][j1] == color).OnlyEnforceIf(edge)
                model.Add(solution[i2][j2] == color).OnlyEnforceIf(edge)
                forward = model.NewBoolVar("")
                backward = model.NewBoolVar("")
                model.AddBoolOr([edge, forward.Not()])
                model.AddBoolOr([edge, backward.Not()])
                model.AddBoolOr([edge.Not(), forward, backward])
                model.AddBoolOr([forward.Not(), backward.Not()])
                arc_variables.append((k1, k2, forward))
                arc_variables.append((k2, k1, backward))
        model.AddCircuit(arc_variables)
    solver = cp_model.CpSolver()
    status = solver.Solve(model)
    if status == cp_model.OPTIMAL:
        for row in solution:
            print("".join(str(solver.Value(x)) for x in row))


if __name__ == "__main__":
    main(
        [
            [1, 0, 0, 2, 3],
            [0, 0, 0, 4, 0],
            [0, 0, 4, 0, 0],
            [0, 2, 3, 0, 5],
            [0, 1, 5, 0, 0],
        ]
    )