Trying to solve a puzzle

OK, I changed loading data from dict to matrix.

I suppose this is the way to declare those variables and constraints:

     @variable(prmx, 
                isCovered[1:_ROWS, 1:_COLUMNS], 
                Bin)

    @variable(prmx, 
                isAntenna[1:_ROWS, 1:_COLUMNS], 
                Bin)

   @objective(
        prmx,
        Min,
        sum(isCovered[i, j] * isHouse[i, j]
            for i in 1:_ROWS, j in 1:_COLUMNS
        )
    )

    @constraint(
        prmx,
        sum(isCovered[i, j] * isHouse[i, j]
            for i in 1:_ROWS, j in 1:_COLUMNS
        )
        >= _TOTAL * 0.7
    )
    

It looks good

Your last set of constraints is useless I think

And the objective should be to minimize the amount of Antennas, you already have the 70% rule as a constraint

Yes, they say “at least 70%” so it has to bi Min with constraint sum(…) >= _TOTAL * 0.7.

Now I can;t quite define constraints for isAntenna → isCovered. Per my logic it should be if there is Antenna it should cover adjusting cells, but still dont know how to write constraints for that.

I think it worked with

    @constraint(
        prmx, 
        [i in 1:_ROWS, j in 1:_COLUMNS],
        isCovered[i, j]
        <=
        isAntenna[i, j] 
        + sum(
            isAntenna[r, c]
            for r in i - 1:i + 1, c in j - 1:j + 1
                if r != i && r >= 1 && r <= _ROWS && c != j && c >= 1 && c <= _COLUMNS
        )   
    )

you can simplify a bit:

    @constraint(
        prmx, 
        [i in 1:_ROWS, j in 1:_COLUMNS],
        isCovered[i, j]
        <=
        sum(
            isAntenna[r, c]
            for r in i - 1:i + 1, c in j - 1:j + 1
            if 1 <= r <= _ROWS && 1 <= c <= _COLUMNS
        )   
    )

and the objective should simply be the total amount of antennas:

@objective(prmx, Min, sum(isAntenna))

then it should be good I think

for @objective I agree and understand, but can’t visualize that isCovered and isAntenna relationship, here it says that isCovered <= sum(isAntenna), hm…, should it be other way around like isCovered >= sum(isAntenna)?

if you have isCovered >= sum(isAntenna), your solver will say,
“ok, isAntenna = 0 everywhere and isCovered = 1 everywhere”

which is of course not want we want (although coverage without antennas would be great :smiley:)

there is proper maths explanations to transform and and or into linear constraints but I couldn’t explain the proper maths theory

I found this post that explains how to transform boolean logic in linear equations

the proper complete way to do a or is (if b is a list):

a = |(b...)

is equivalent to

a <= sum(b)
a .>= b

in this case the a .>= b isn’t strictly needed (I think) since the solver will try to minimize b and maximize a anyway

What do you think about this:

    @constraint(
        prmx, 
        [i in 2:_ROWS-1, j in 2:_COLUMNS-1],
          sum(
            isCovered[r, c]
            for r in i - 1:i + 1, c in j - 1:j + 1
          )  - 9 * isAntenna[i, j]  >= 0
    )
isCovered[i, j] <=
sum(
    isAntenna[r, c]
    for r in i-1:i+1, c in j-1:j+1
    if 1 <= r <= _ROWS && 1 <= c <= _COLUMNS
) 

corresponds to

isCovered[i, j] ==
|([
    isAntenna[r, c]
    for r in i-1:i+1, c in j-1:j+1
    if 1 <= r <= _ROWS && 1 <= c <= _COLUMNS
]...) 

and for completeness you can add:

isCovered[i, j] .>=
[
    isAntenna[r, c]
    for r in i-1:i+1, c in j-1:j+1
    if 1 <= r <= _ROWS && 1 <= c <= _COLUMNS
]

this works fine for me and I find the right answer (6 antennas for 29 houses covered)
I don’t understand what you’re trying to do in your last post

1 Like

Here’s a more compact way to write it:

using JuMP, Cbc
function main(; houses::Vector{Tuple{Int,Int}}, ratio::Float64)
    M = maximum(h[1] for h in houses)
    N = maximum(h[2] for h in houses)
    model = Model(Cbc.Optimizer)
    @variable(model, cell_tower[1:M, 1:N], Bin)
    @variable(model, is_covered[houses], Bin)
    @objective(model, Min, sum(cell_tower))
    stencil(i, j) = cell_tower[max(i-1,1):min(i+1,M), max(j-1,1):min(j+1,N)]
    @constraint(model, [h in houses], is_covered[h] <= sum(stencil(h[1], h[2])))
    @constraint(model, sum(is_covered) >= ratio * length(houses))
    optimize!(model)
    return value.(cell_tower)
end
houses = [(1, 9), (3, 2), (4, 4), (5, 6), (8, 9), (9, 1)]
main(houses, ratio = 0.7)
2 Likes

The most important, what I failed to understand, is that isCovered is in the function of houses. I was thinking about the grid that one antenna can cover, but I should think about which houses antenna(s) can cover instead.

Here is the model:

using JuMP, Cbc, XLSX

mat  = XLSX.readdata("data/input.xlsx", "sheet1", "A1:J10")

M = size(mat, 1)
N = size(mat, 2)
ratio = 0.7

houses = []
for i = 1 : M
  for j = 1 : N
    if mat[i,j] == 1
      push!(houses, (i,j))
    end
  end
end

model = Model(Cbc.Optimizer)

@variable(model, 
            isCovered[houses], 
            Bin)

@variable(model, 
            isAntenna[1:M, 1:N], 
            Bin)

@objective(
    model,
    Min,
    sum(isAntenna)
)

@constraint(
    model,
    sum(isCovered) >= length(houses) * ratio
)

@constraint(
    model, 
    [h in houses],
    isCovered[h] <=
    sum(
        isAntenna[r, c]
        for r in h[1]-1:h[1]+1, c in h[2]-1:h[2]+1
        if 1 <= r <= M && 1 <= c <= N
    ) 
)

optimize!(model)

res = primal_status(model)

if res == MOI.FEASIBLE_POINT
    val = value.(isAntenna)
    open("out/result.txt","w") do io
        for i in 1:M
            for j in 1:N
                print(io, floor(Int, val[i,j]))
                print(io, '\t')
            end
            print(io, "\n")
        end
    end
end

Got it finally! 6 antennas → 29 houses

Great!

To see a Matrix approach, here was my take on the problem:

using JuMP, Cbc
m = Model(Cbc.Optimizer)
set_silent(m)

houses = [
0 0 1 0 1 1 0 1 0 0
1 0 0 0 1 0 0 0 0 0
0 0 0 1 0 1 1 0 1 1
0 0 0 0 1 0 1 0 1 0
0 0 1 0 0 1 1 1 1 0
1 0 0 1 1 0 1 0 0 0
1 1 1 1 1 0 0 0 0 1
0 0 0 0 1 0 1 0 0 0
1 0 1 0 0 0 1 1 0 1
1 0 0 0 1 1 0 0 1 1
]

@variables m begin 
    isantenna[1:12,1:12], Bin
    iscovered[1:10,1:10], Bin
end

@constraints m begin 
    isantenna[  1,  :] .== 0
    isantenna[end,  :] .== 0
    isantenna[  :,  1] .== 0
    isantenna[  :,end] .== 0
        
    iscovered .<= sum(isantenna[i:end-3+i, j:end-3+j] for i in 1:3, j in 1:3)

    sum(iscovered .* houses) >= 0.7*sum(houses)
end

@objective m Min sum(isantenna)

optimize!(m)

result = Int.(value.(isantenna[2:end-1, 2:end-1]))
println(sum(result))
display(result)
1 Like

Great! But I am still a conventional programmer though :smile: and will need some time to understand and start using those operators like . Anyway I am already on another page: puzzlor.com did you solve that too?

hmmm, that’s not linear (and possibly not convex)

I might try it when I have some time

1 Like

with 1,000,000 possibilities only, it’s totaly reasonable to brute-force the problem
without optimizing my code, I find the solution in 7 seconds

2 Likes

Wait, you already created a model for this problem and solve it???

not a JuMP model, you can’t express this as a linear problem, and I don’t think it’s even a convex problem

just a few loops over the 1,000,000 possibilities and keeping the best arrangement

1 Like

I was thinking is it possible to put a @variable in for loop? For example:


Size = 2

@constraint(
    model, 
    [h in houses],
    isCovered[h] <=
    sum(
        isAntenna[r, c]
        for r in h[1] - Size:h[1] + Size, c in h[2] - Size:h[2] + Size
        if 1 <= r <= M && 1 <= c <= N
    ) 
)

Is it possible to make Size to be Int @variable? ...and if not, why?

Is it possible to make Size to be Int @variable?

No.

if not, why?

Sums have to be structural over a fixed number of elements. The values of the variables can change, but not the sparsity pattern of the summation.