Here’s how I would write it, using the package GridGraphs.jl that I developed (be sure to install the latest version 0.9.1):

```
using Graphs
using GridGraphs: GridGraph, QUEEN_DIRECTIONS, index_to_coord
using HiGHS
using JuMP
d = 4 # size of the grid, vertices are named u or v
n = 9 # maximum number, numbers are named k
grid = GridGraph(ones(Int, d, d); directions=QUEEN_DIRECTIONS);
numbers = rand(1:n, d^2); # sequence of numbers (order unimportant)
numbers_count = [count(==(k), numbers) for k in 1:n];
model = Model(HiGHS.Optimizer)
# Variables: x[v, k] = 1 <-> number k on vertex v
@variable(model, x[v=1:d^2, k=1:n], Bin);
# Auxiliary variables: y[u, v, k] = x[u, k] x[v, k]
@variable(model, y[u=1:d^2, v=1:d^2, k=1:n; has_edge(grid, u, v)], Bin);
# Constraint 1: one number per vertex
@constraint(model, [v = 1:d^2], sum(x[v, :]) == 1);
# Constraint 2: as many k's as needed
@constraint(model, [k = 1:n], sum(x[:, k]) == numbers_count[k]);
# Linearization constraints: coherence between y and x
@constraint(model, [u = 1:d^2, v = 1:d^2, k = 1:n; has_edge(grid, u, v)], y[u, v, k] <= x[u, k]);
@constraint(model, [u = 1:d^2, v = 1:d^2, k = 1:n; has_edge(grid, u, v)], y[u, v, k] <= x[v, k]);
@constraint(model, [u = 1:d^2, v = 1:d^2, k = 1:n; has_edge(grid, u, v)], y[u, v, k] >= x[u, k] + x[v, k] - 1);
# Objective: sum of k xᵤₖ xᵥₖ over all k's and all edges (u, v)
@objective(model, Max, sum(k * y[src(e), dst(e), k] for k in 1:n, e in edges(grid)));
optimize!(model)
x_sol = round.(Int, value.(x))
solution = zeros(Int, d, d)
for v in vertices(grid)
(i, j) = index_to_coord(grid, v)
solution[i, j] = argmax(x_sol[v, :])
end
```

Here’s what we get at the end for one random instance:

```
julia> 1:n .=> numbers_count
9-element Vector{Pair{Int64, Int64}}:
1 => 0
2 => 2
3 => 1
4 => 2
5 => 4
6 => 0
7 => 3
8 => 3
9 => 1
julia> solution
4×4 Matrix{Int64}:
4 4 3 9
8 8 5 5
7 8 5 5
7 7 2 2
```

It takes a VERY long time to solve for `d=5`

so maybe there is a better formulation though