Okay, I tried your suggested code with Gurobi and get the following output:
cademic license - for non-commercial use only - expires 2022-10-27
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 242 rows, 484 columns and 484 nonzeros
Model fingerprint: 0x2b837d1e
Model has 121 quadratic objective terms
Model has 121 quadratic constraints
Variable types: 121 continuous, 363 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
QMatrix range [1e+00, 1e+00]
QLMatrix range [1e+00, 8e+01]
Objective range [0e+00, 0e+00]
QObjective range [2e+00, 2e+00]
Bounds range [2e+01, 2e+01]
RHS range [2e+01, 2e+01]
QRHS range [3e+01, 2e+03]
Presolve removed 138 rows and 276 columns
Presolve time: 0.30s
Presolved: 260 rows, 208 columns, 676 nonzeros
Presolved model has 52 quadratic objective terms
Presolved model has 52 bilinear constraint(s)
Variable types: 52 continuous, 156 integer (0 binary)
Root relaxation: objective 8.917000e+03, 451 iterations, 0.00 seconds
H 3676 3984 11745.000000 8917.00000 24.1% 2.8 1s
H 4017 3019 10557.000000 8917.00000 15.5% 2.8 1s
* 4024 2851 798 10494.000000 8917.00000 15.0% 2.8 1s
H 4601 2389 9864.0000000 8917.00000 9.60% 2.7 1s
H 4601 2269 9801.0000000 9721.00000 0.82% 2.7 1s
Cutting planes:
MIR: 1
Explored 4641 nodes (13520 simplex iterations) in 1.71 seconds
Thread count was 8 (of 8 available processors)
Solution count 6: 9801 9864 10494 ... 11808
Optimal solution found (tolerance 1.00e-04)
Best objective 9.801000000000e+03, best bound 9.801000000000e+03, gap 0.0000%
User-callback calls 9456, time in user-callback 0.00 sec
So far so good. But when I look at the output for a simple test case, it does not appear well-converged.
Left is z′, middle is manually wrapping the values to within 0:16, and right is z after optimization:

Are there any options for Gurobi that I should tune for this type of problem?
Thanks very much!