The urban planning problem

An "urban planning" problem based on an example from puzzlor.

using JuMP
import HiGHS
import Test

function example_urban_plan()
    model = Model(HiGHS.Optimizer)
    # x is indexed by row and column
    @variable(model, 0 <= x[1:5, 1:5] <= 1, Int)
    # y is indexed by R or C, the points, and an index in 1:5. Note how JuMP
    # allows indexing on arbitrary sets.
    rowcol = ["R", "C"]
    points = [5, 4, 3, -3, -4, -5]
    @variable(model, 0 <= y[rowcol, points, 1:5] <= 1, Int)
    # Objective - combine the positive and negative parts
    @objective(
        model,
        Max,
        sum(
            3 * (y["R", 3, i] + y["C", 3, i]) +
            1 * (y["R", 4, i] + y["C", 4, i]) +
            1 * (y["R", 5, i] + y["C", 5, i]) -
            3 * (y["R", -3, i] + y["C", -3, i]) -
            1 * (y["R", -4, i] + y["C", -4, i]) -
            1 * (y["R", -5, i] + y["C", -5, i]) for i in 1:5
        )
    )
    # Constrain the number of residential lots
    @constraint(model, sum(x) == 12)
    # Add the constraints that link the auxiliary y variables to the x variables
    for i in 1:5
        @constraints(model, begin
            # Rows
            y["R", 5, i] <= 1 / 5 * sum(x[i, :]) # sum = 5
            y["R", 4, i] <= 1 / 4 * sum(x[i, :]) # sum = 4
            y["R", 3, i] <= 1 / 3 * sum(x[i, :]) # sum = 3
            y["R", -3, i] >= 1 - 1 / 3 * sum(x[i, :]) # sum = 2
            y["R", -4, i] >= 1 - 1 / 2 * sum(x[i, :]) # sum = 1
            y["R", -5, i] >= 1 - 1 / 1 * sum(x[i, :]) # sum = 0
            # Columns
            y["C", 5, i] <= 1 / 5 * sum(x[:, i]) # sum = 5
            y["C", 4, i] <= 1 / 4 * sum(x[:, i]) # sum = 4
            y["C", 3, i] <= 1 / 3 * sum(x[:, i]) # sum = 3
            y["C", -3, i] >= 1 - 1 / 3 * sum(x[:, i]) # sum = 2
            y["C", -4, i] >= 1 - 1 / 2 * sum(x[:, i]) # sum = 1
            y["C", -5, i] >= 1 - 1 / 1 * sum(x[:, i]) # sum = 0
        end)
    end
    # Solve it
    optimize!(model)
    Test.@test termination_status(model) == OPTIMAL
    Test.@test primal_status(model) == FEASIBLE_POINT
    Test.@test objective_value(model) ≈ 14.0
    return
end

example_urban_plan()

Presolving model
61 rows, 85 cols, 385 nonzeros
61 rows, 85 cols, 385 nonzeros
Objective function is integral with scale 1

Solving MIP model with:
   61 rows
   85 cols (85 binary, 0 integer, 0 implied int., 0 continuous)
   385 nonzeros

( 0.0s) Starting symmetry detection
( 0.0s) Found 5 generators

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   50              -inf                 inf        0      0      0         0     0.0s
 S       0       0         0   0.00%   50              -6               933.33%        0      0      0         0     0.0s
         0       0         0   0.00%   28.8            -6               580.00%        0      0      4        63     0.0s
 H       0       0         0   0.00%   14              14                 0.00%     3645     99     34       642     0.2s

Solving report
  Status            Optimal
  Primal bound      14
  Dual bound        14
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    14 (objective)
                    0 (bound viol.)
                    2.22044604925e-16 (int. viol.)
                    0 (row viol.)
  Timing            0.16 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             1
  LP iterations     642 (total)
                    0 (strong br.)
                    579 (separation)
                    0 (heuristics)