What strategy to allocate resources while respecting many constraints

rocco_sprmnt21 · May 27, 2023, 2:16pm

Suppose we have a table with 8 rows and 30 (or 31) columns.
We have to label each cell with G and N, respecting the following rules (sometimes some other constraint is added )

Each column (i.e. each day) must have exactly one G and one N, eventually in the same cell (GN)
Each row must have exactly g G and n N, where g and n are given input values, usually in {3,4} [of course sum(g)==sum(n)==30(31)]
The cells containing xx are forbidden for both N and G; xg are forbidden for G and xn are forbidden for N
For each row, 3 consecutive days with G (or with N) for the same row are forbidden
For each row, 3 days with GN are forbidden (preferably only one GN, except for row r, “who” prefer as many GNs as possible)
For each row you can’t have G following a N
Can’t have 3 cells ( N or G) during the weeks end in the same row

PS
I solved the practical problem in an “artisanal” way, but I’m curious to know if this falls into some category that is treated by some package or in any case algorithms suitable for this type of problem are known

I attach an example of an input table

	gio	ven	sab	dom	lun	mar	mer	gio	ven	sab	dom	lun	mar	mer	gio	ven	sab	dom	lun	mar	mer	gio	ven	sab	dom	lun	mar	mer	gio	ven					
	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30				g	n
R1	xx	xx	xx	xx					xx								xx																	4	3
R2	xx	xx	xx	xx	xx	xx	xx								xg	xx	xx	xx	xx										xx	xx				4	4
R3		xx		xg	xx	xx	xx	xx	xx	xx	xx		xx	xx	xx	xx	xx	xg			xx	xx	xx	xx				xx	xx	xx				4	4
R4	xg			xg		xx	xx			xx	xx						xx	xx						xx	xx									4	3
R5		xx	xx	xx			xx		xg	xg						xg	xx	xx	xx	xx	xx	xx	xx	xx	xx	xx	xx	xx	xx	xx				3	4
R6	xx	xx	xx	xx	xn	xg	xn		xx			xx		xn		xn			xx		xn		xx	xx	xx	xn		xn		xn				4	4
R7	xx	xx	xx																				xx	xx	xx	xx	xx	xx	xx	xx				3	4
R8	xx	xx	xx	xx	xx											xx	xx	xx					xx	xx	xx			xx	xx	xx				4	4

DrChainsaw · May 27, 2023, 4:43pm

I haven’t read the constraints very carefully, but it sounds like something JuMP would allow you to do easily. I would use one 30x8 set of binary variables for G and one for N and let 1 mean that there is a G/N in that position. Check out the soduko example for a start.

rocco_sprmnt21 · May 28, 2023, 6:28am

Thank you.
It seems that this is a package that can be used for “my” problem.
To practice a bit I tried to apply it to the classic case of the 8 queens, which (to my eye) seems more similar to “my” problem.

using JuMP
using HiGHS

queen = Model(HiGHS.Optimizer)
set_silent(queen)
@variable(queen, x[i = 1:8, j = 1:8], Bin)


for ind in 1:8  
        @constraint(queen, sum(x[ind, j] for j in 1:8) == 1)
        @constraint(queen, sum(x[i, ind] for i in 1:8) == 1)
end
for ind in 2:16
    @constraint(queen,  sum(x[r, ind-r]  for r in 1:ind-1 if 0<r<9 && 0<ind-r<9) <= 1)
end

for ind in 7:-1:-7
@constraint(queen,  sum(x[ind+r,r]  for r in 1:8-ind if 0<r<9 && 0<ind+r<9) <= 1)
end


optimize!(queen)
csp_sol = round.(Int, value.(x))

PS
in this case what is the optimized (minimized?) objective function?

rocco_sprmnt21 · May 28, 2023, 9:06am

I have a(n unsatisfactory) sol of this type, where in the row=1 there are 2 GN

julia> sol
8×30 Matrix{String}:
 ""   ""   "G"  ""   "X"  ""   "X"  "N"  ""   "X"  "X"   "X"  ""   "GN"  "X"  ""   ""    "X"  "X"  ""   "GN"  "X"   "X"  "X"   "G"  "N"  "X"  "X"  "X"  "X"
 "X"  "X"  "N"  "X"  "X"  "G"  "N"  "X"  "G"  "X"  "X"   ""   "X"  "X"   ""   "X"  "GN"  "X"  ""   "X"  "X"   ""    ""   "X"   "N"  "G"  "X"  "X"  "X"  ""
 ""   "G"  ""   "X"  "X"  ""   "X"  "X"  "X"  "G"  "X"   "X"  ""   ""    "X"  "N"  "X"   "N"  ""   "G"  "X"   ""    "N"  "X"   "X"  ""   ""   "N"  ""   "X"
 "G"  ""   "X"  "X"  "X"  "X"  "X"  "G"  "X"  "N"  "GN"  "N"  ""   "X"   "X"  ""   "X"   ""   ""   "N"  "X"   "X"   ""   "X"   "X"  "X"  "X"  "X"  ""   "X"
 "X"  ""   "X"  "X"  "N"  "N"  "G"  ""   "X"  "X"  ""    "X"  "G"  ""    "X"  "X"  "X"   ""   "X"  "X"  ""    "GN"  ""   "X"   "X"  "X"  "X"  ""   "G"  "X"
 ""   "X"  "X"  "N"  "X"  "X"  "X"  ""   "X"  "X"  ""    "G"  "N"  "X"   "X"  "X"  ""    "G"  "G"  "X"  ""    "X"   "G"  "X"   "X"  "X"  ""   ""   "N"  "X"
 ""   "N"  "X"  "G"  "X"  "X"  ""   ""   "N"  ""   "X"   "X"  ""   ""    "G"  "G"  ""    "X"  "N"  ""   "X"   "X"   ""   "X"   "X"  ""   "N"  "X"  "X"  "G"
 "N"  "X"  ""   ""   "G"  "X"  ""   "X"  "X"  "X"  ""    "X"  ""   "X"   "N"  ""   "X"   "X"  "X"  "X"  "X"   ""    "X"  "GN"  ""   "X"  "G"  "G"  ""   "N"

I tried, but without success, to enforce such a constraint

for r in 1:8
    @constraint(indisp, sum(x[r, g, 1]*x[r, g, 2] for g in 1:30) <=2 )
end
#ERROR: Constraints of type MathOptInterface.ScalarQuadraticFunction{Float64}-in-MathOptInterface.LessThan{Float64} are not supported by the solver.

# also this form get an error

julia> for r in 1:8
           @constraint(indisp, sum(x[r, g, 1]&&x[r, g, 2] for g in 1:30) <=2 )
       end
ERROR: TypeError: non-boolean (VariableRef) used in boolean context

the initial part of the code

indisp = Model(HiGHS.Optimizer)
set_silent(indisp)
@variable(indisp, x[i = 1:8, j = 1:30, k=1:2], Bin)
M=[4 4; 4 4; 3 4; 3 4; 4 3; 4 3; 4 4; 4 4]

# turni nel mese per operatore
for k in 1:2  
    for r in 1:8  
        @constraint(indisp, sum(x[r, j, k] for j in 1:30) == M[r,k])
    end
end
# turni nel giorno 

for k in 1:2  
    for g in 1:30
        @constraint(indisp, sum(x[r, g, k] for r in 1:8) == 1)
    end
end

rocco_sprmnt21 · May 28, 2023, 10:43am

using JuMP
using HiGHS

indisp = Model(HiGHS.Optimizer)
set_silent(indisp)
@variable(indisp, x[i = 1:8, j = 1:30, k=1:2], Bin)
M=[4 4; 4 4; 3 4; 3 4; 4 3; 4 3; 4 4; 4 4]

# turni nel mese per operatore
for k in 1:2  
    for r in 1:8  
        @constraint(indisp, sum(x[r, j, k] for j in 1:30) == M[r,k])
    end
end

# each day must have just one k==1 and one k==2

for k in 1:2  
    for g in 1:30
        @constraint(indisp, sum(x[r, g, k] for r in 1:8) == 1)
    end
end

# for each row  max a day having  both  G and N
# for r in 1:8
#     @constraint(indisp, sum((x[r, g, 1] && x[r, g, 2]) for g in 1:30) <=2 )
# end
# tried to simulate the above constraint in the following way

for r in 1:8
    @constraint(indisp, sum(x[r, g, 1] + x[r, g, 2]-1 for g in 1:30) <=2 )
end

init_sol = [rand(Bool) for r in 1:8, g in 1:30]

for i in 1:8  # fix to 0 the 'forbidden' cells
    for j in 1:30
        if init_sol[i, j] == 0
            for k in 1:2
            fix(x[i, j, k], 0; force = true)
            end
        end
    end
end

optimize!(indisp)

x_val = value.(x)

sol = ["" for r in 1:8, g in 1:30]

for r in 1:8
    for g in 1:30
        if !init_sol[r,g]
            sol[r,g]="X"
        else
            sol[r,g]=""
        end
    end
end

for i in 1:8
    for j in 1:30
        for k in 1:2
            if round(Int, x_val[i, j, k]) == 1
                sol[i, j] *= k==1 ? 'G' : 'N'
            end
        end
    end
end
sol

But … this trick doesn’t work either

julia> sol
8×30 Matrix{String}:
 "X"   ""    ""    "G"  "X"   "X"   "X"  "N"  ""    "G"  ""    ""   ""    "X"  "G"  …  ""    ""    "N"  "X"   "X"  "N"  "X"  "G"  ""   ""    "X"   "X"  "N"  "X"     
 "X"   "X"   ""    ""   ""    ""    "N"  ""   "X"   "X"  "X"   "G"  "GN"  "X"  "X"     "X"   ""    "X"  "X"   "X"  "X"  ""   "N"  "X"  "X"   "X"   ""   "G"  "GN"    
 "X"   "X"   "GN"  "X"  "X"   "X"   "X"  ""   "X"   "X"  "X"   "N"  "X"   "N"  "X"     "X"   "GN"  "X"  "X"   "X"  "X"  "G"  ""   "X"  "X"   "X"   "X"  ""   ""      
 "X"   "X"   "X"   "N"  "X"   "X"   ""   "X"  "X"   ""   "X"   ""   ""    ""   ""      "X"   "X"   "X"  ""    "G"  "X"  "N"  "X"  "N"  "GN"  "X"   "G"  ""   "X"     
 "X"   "X"   "X"   "X"  "X"   "X"   "G"  "X"  "X"   "X"  "X"   "X"  "X"   "G"  "N"     ""    ""    ""   "X"   "N"  "G"  "X"  "X"  "G"  "X"   ""    "X"  "X"  ""
 "GN"  ""    "X"   "X"  "X"   "X"   ""   "G"  "X"   "X"  ""    ""   "X"   "X"  "X"  …  ""    ""    "G"  "GN"  ""   "X"  "X"  "X"  "X"  "X"   "X"   "N"  "X"  "X"     
 "X"   "X"   ""    ""   "GN"  "GN"  ""   "X"  "GN"  ""   ""    ""   ""    ""   ""      "X"   "X"   "X"  ""    ""   "X"  "X"  "X"  "X"  "X"   "GN"  "X"  "X"  "X"     
 "X"   "GN"  ""    "X"  "X"   ""    "X"  "X"  "X"   "N"  "GN"  "X"  "X"   "X"  "X"     "GN"  "X"   "X"  ""    "X"  "X"  "X"  "X"  "X"  "X"   "X"   "X"  "X"  ""

DrChainsaw · May 28, 2023, 1:03pm

Gettng mixed integer progamming to do what you want can certainly be a bit of a mind-warping puzzle, and somethings just don’t work. A good high level tip is that it is searchable, often by by search terms like “ integer programming”.

For the particular problem, just as you noted, you want boolean and between two variables so you can put a constraint on that. This has a linear formulation by introducing a new variable and putting some clever constraints on it.

The following (untested) function returns a new variable which is the element wise and of var1 and var2.

function andvar!(model, var1, var2)
    andvar = @variable(model, [1:length(var1)], Bin)
    @constraint(model, andvar .>= var1 .+ var2 .- 1)
    @constraint(model, andvar .<= var1)
    @constraint(model, andvar .<= var2)
    return andvar
end

you can then constrain the sum of this variable to be less than or equal to 3 for each row.

I think you might run into similar issues with constraint 4 as the way to write it might not be intuitive, but it is a common enough problem that it should show up in a search (e.g. something like “contiguous binary variables integer programming”). Here is a function which does something similar, but it does a couple of other things too, so it might be easier to just write it from scratch.

As for selecting an objective. You might not need one, but one idea is to try to formulate things like “preferably only one GN” in the objective, perhaps just by minimizing the variable which counts the number of GNs that you anyways have created for the constraint.

rocco_sprmnt21 · May 28, 2023, 1:48pm

Thank you.
I will study your answer carefully.
In the meantime I’ve found a way to use some “trick” suggested by the manual to use AND and OR
for rules 4. and 5.

using JuMP
using HiGHS

indisp = Model(HiGHS.Optimizer)
set_silent(indisp)

# k==1 is for G, k==2 is for N, k==3 is for G/\N, k==4 is for G\/N

@variable(indisp, x[i = 1:8, j = 1:30, k=1:4], Bin)
M=[4 4; 4 4; 3 4; 3 4; 4 3; 4 3; 4 4; 4  4]

# totsl of G and N for row
for k in 1:2  
    for r in 1:8  
        @constraint(indisp, sum(x[r, j, k] for j in 1:30) == M[r,k])
    end
end

# each day maust have exactly one G (k==1) and one N (k==2)
for k in 1:2  
    for g in 1:30
        @constraint(indisp, sum(x[r, g, k] for r in 1:8) == 1)
    end
end


#  x[3] = x[1]/\x[2]
for r in 1:8, j in 1:30
    @constraints(indisp, begin
       x[r,j,3] <= x[r,j,1]
       x[r,j,3] <= x[r,j,2]           
       x[r,j,3] >= x[r,j,1] + x[r,j,2] - 1
    end)
end
#max  one GN for row
for r in 1:8
    @constraint(indisp, sum(x[r, g, 3]  for g in 1:30) <= 1 )
end

# no more than 2 consecutive G[N], G[N]  
#  x[4] = x[1]\/x[2]
for r in 1:8, j in 1:30 
    @constraints(indisp, begin
           x[r,j,1] <= x[r,j,4]
           x[r,j,2] <= x[r,j,4]
           x[r,j,4] <= x[r,j,1] + x[r,j,2]
    end)
end

for r in 1:8
    for g in 3:30
         @constraint(indisp, sum(x[r, i, 4] for i in g-2:g) <= 2 )
    end
end

init_sol = [rand(Bool) for r in 1:8, g in 1:30] 

for i in 1:8
    for j in 1:30
        if init_sol[i, j] == 0
            for k in 1:2
            fix(x[i, j, k], 0; force = true)
            end
        end
    end
end



optimize!(indisp)

x_val = value.(x)

sol = ["" for r in 1:8, g in 1:30]
for r in 1:8
    for g in 1:30
        if !init_sol[r,g]
            sol[r,g]="X"
        else
            sol[r,g]=""
        end
    end
end

for i in 1:8
    for j in 1:30
        for k in 1:2
            if round(Int, x_val[i, j, k]) == 1
                sol[i, j] *= k==1 ? 'G' : 'N'
            end
        end
    end
end
sol

julia> sol
8×30 Matrix{String}:
 "N"  "X"  "X"  "X"   ""   ""   "X"  "X"  "X"  "X"   "X"  "G"  "G"  ""    "X"  "X"  "N"  "X"  "N"  "X"  "GN"  ""    "X"  "X"  ""    "X"  "X"  "G"  ""   "X"
 "X"  "G"  ""   "X"   "N"  ""   "X"  "N"  "G"  "X"   "X"  ""   "X"  ""    "G"  "X"  "X"  "N"  "X"  "X"  ""    "X"   "X"  ""   "X"   ""   ""   "X"  "G"  "N"
 ""   ""   ""   "GN"  "X"  ""   ""   ""   "N"  "X"   "G"  ""   ""   "X"   "X"  "G"  "X"  "X"  ""   "X"  "X"   ""    ""   "X"  "X"   "N"  ""   ""   "N"  "X"
 "X"  "X"  "X"  "X"   "X"  ""   "N"  "X"  "X"  "X"   ""   ""   "X"  "X"   "X"  "N"  "X"  ""   "G"  "X"  "X"   ""    ""   "G"  "GN"  "X"  "N"  "X"  ""   "X"
 "X"  "N"  "X"  ""    ""   "X"  "G"  "X"  "X"  "X"   ""   "N"  "X"  "GN"  "X"  "X"  "X"  "X"  "X"  "X"  "X"   ""    ""   "X"  ""    "X"  "G"  "X"  "X"  "G"
 ""   "X"  "G"  ""    "G"  ""   "X"  "X"  "X"  "GN"  "X"  ""   "X"  "X"   "N"  "X"  "G"  "X"  "X"  ""   "X"   ""    ""   "N"  "X"   ""   ""   "X"  ""   "X"
 "G"  "X"  "N"  ""    "X"  "G"  ""   "X"  "X"  "X"   "N"  "X"  "N"  ""    "X"  "X"  "X"  "X"  "X"  "G"  "X"   "X"   "G"  "X"  ""    ""   "X"  "N"  ""   "X"
 "X"  ""   "X"  "X"   ""   "N"  ""   "G"  "X"  "X"   ""   ""   "X"  ""    ""   "X"  ""   "G"  "X"  "N"  ""    "GN"  "N"  ""   ""    "G"  ""   "X"  "X"  ""

rocco_sprmnt21 · May 28, 2023, 6:21pm

Thanks to @DrChainsaw .@odow

This should be the solution that satisfies all the rules, expected.
At a later stage I will add some comments (or some slight modifications) to make clearer (starting with myself, in some time) the meaning of the various expressions.

(edited after odow’s suggestions)

using JuMP
import HiGHS
M = Dict("G" => [4, 4, 3, 3, 4, 4, 4, 4], "N" => [4, 4, 4, 4, 3, 3, 4, 4])
startday = 7
we = [6, 7, 13, 14, 20, 21, 27, 28, 34, 35]
wse = filter(d-> 0 < d < 32, we .- (startday - 1))
init_sol = [rand(Bool) for r in 1:8, g in 1:30] 
init_sol_1 = [(rand(1:8), rand(1:30)) for _ in 1:10]
init_sol_2 = [(rand(1:8), rand(1:30)) for _ in 1:10]
model = Model(HiGHS.Optimizer)
# set_silent(model)
@variable(model, x[i in 1:8, j in 1:30, k in ["G", "N", "G/\\N", "G\\/N"]], Bin)
@constraints(model, begin
    # total of G and N for row
    [k in ["G", "N"], r in 1:8], sum(x[r,:,k]) == M[k][r]
    # each day must have exactly one G (k==1) and one N (k==2)
    [k in ["G", "N"], g in 1:30], sum(x[:,g,k]) == 1
    # max of G and N for row for weekends
    #[k in ["G", "N", "G/\\N"], r in 1:8], sum(x[r,j,k] for j in wse) <= 3
    [ r in 1:8], sum(x[r,j,"G\\/N"] for j in wse) <= 3
    # max one GN for row
    [r in 1:8], sum(x[r,g,"G/\\N"] for g in 1:30) <= 1
    #  x["G/\N"] = x["G"] /\ x["N"]
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] <= x[r,j,"G"]
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] <= x[r,j,"N"]           
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] >= x[r,j,"G"] + x[r,j,"N"] - 1

    #  x["G\/N"] = x["G"] \/ x["N"]\/x["G/\\N"]
    [r in 1:8, j in 1:30], x[r,j,"G"] <= x[r,j,"G\\/N"]
    [r in 1:8, j in 1:30], x[r,j,"N"] <= x[r,j,"G\\/N"]
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] <= x[r,j,"G\\/N"]
    [r in 1:8, j in 1:30], x[r,j,"G\\/N"] <= x[r,j,"G"] + x[r,j,"N"]+ x[r,j,"G/\\N"]
    # no more than 2 elements consecutive from {G, N, GN} 
    [r in 1:8, g in 3:30], sum(x[r,i,"G\\/N"] for i in g-2:g) <= 2
    # no G can follow N 
    [r in 1:8, g in 2:30], x[r,g,"G"] <= 1 - x[r,g-1,"N"]
    [r in 1:8, g in 2:30], x[r,g,"G"] <= 1 - x[r,g-1,"G/\\N"]
    [r in 1:8, g in 2:30], x[r,g,"G/\\N"] <= 1 - x[r,g-1,"N"]
    [r in 1:8, g in 2:30], x[r,g,"G/\\N"] <= 1 - x[r,g-1,"G/\\N"]
end)
sol = [init_sol[r, g] ? "." : "X" for r in 1:8, g in 1:30]
for i in 1:8, j in 1:30
    if init_sol[i, j] == 0
        fix.(x[i, j, ["G", "N"]], 0; force = true)
    end
end
for (i, j) in init_sol_1
    fix(x[i, j, "G"], 0; force = true)
    sol[i, j] = "xg"
end
for (i, j) in init_sol_2
    fix(x[i, j, "N"], 0; force = true)
    sol[i, j] = "xn"
end
optimize!(model)
for i in 1:8, j in 1:30, k in ["G", "N"]
    if round(Int, value(x[i, j, k])) == 1
        sol[i, j] *= k
    end
end
for i in 1:8
    for j in 1:30
        print(lpad(sol[i, j], 4))
    end
    println()
end

.G   X   .   X  .N   .   . .GN   .   X   .   .   X   X   X  .N   X   X   X xgN   X   .   X   .  .G  .G   .   .   .   .
  xn   . xnG   X  xg  xn   X  xg  .G   .   X   X   X   X .GN  xn   X  .N   X   .   X   .   X   X   X  .N   X  .G xgN   X
   .   X  .N  xg  .G  .N   .   X  xn   .   X   .  xg   X   X   .   X   X   X   X .GN   .  .N   X   .   .   X   X  .G   X
   .  xg  xg   .   .   . .GN   X   X   X  xn  .N  .N   X   X   X   X  xg   X   X   .  .G   .   X   X   X  .G  .N  xn   X
   .   X   X   .   X   X   X   X   X  .G   X  .G   .  .G   X   X .GN   X  .N   .   .  .N   X   X   .   X   X   X   X   .
  .N   .   X  .G   .  .G   .   X  .N   .   X   X  .G   X   .   .   X  .G   X   .   X   .   .   X   .   X  .N   .   X   X
   X  .N   X  .N   .   X   .   X   X  .N  .N   X   X   X   X  .G   X   X   X  .G   X   X  .G   X   X   X   X   X   X  .G
   .  .G   X   .   X   .  xn   X   .   X  .G   .   X  .N   .  xn   X   X xnG   X   .   .   X .GN xgN   .   X   X   .  .N

odow · May 28, 2023, 7:56pm

One thing that would make it clearer is using strings for the k indices. Here’s another way of writing your model:

using JuMP
import HiGHS
M = Dict("G" => [4, 4, 3, 3, 4, 4, 4, 4], "N" => [4, 4, 4, 4, 3, 3, 4, 4])
sd = 7
we = [6, 7, 13, 14, 20, 21, 27, 28, 34, 35]
wse = filter(d-> 0 < d < 32, we .- (sd - 1))
init_sol = [rand(Bool) for r in 1:8, g in 1:30] 
init_sol_1 = [(rand(1:8), rand(1:30)) for _ in 1:10]
init_sol_2 = [(rand(1:8), rand(1:30)) for _ in 1:10]
model = Model(HiGHS.Optimizer)
# set_silent(model)
@variable(model, x[i in 1:8, j in 1:30, k in ["G", "N", "G/\\N", "G\\/N"]], Bin)
@constraints(model, begin
    # total of G and N for row
    [k in ["G", "N"], r in 1:8], sum(x[r,:,k]) == M[k][r]
    # each day must have exactly one G (k==1) and one N (k==2)
    [k in ["G", "N"], g in 1:30], sum(x[:,g,k]) == 1
    # max of G and N for row for weeks end
    [k in ["G", "N"], r in 1:8], sum(x[r,j,k] for j in wse) <= 3
    # max one GN for row
    [r in 1:8], sum(x[r,g,"G/\\N"] for g in 1:30) <= 1
    #  x["G/\N"] = x["G"] /\ x["N"]
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] <= x[r,j,"G"]
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] <= x[r,j,"N"]           
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] >= x[r,j,"G"] + x[r,j,"N"] - 1
    #  x["G\/N"] = x["G"] \/ x["N"]
    [r in 1:8, j in 1:30], x[r,j,"G"] <= x[r,j,"G\\/N"]
    [r in 1:8, j in 1:30], x[r,j,"N"] <= x[r,j,"G\\/N"]
    [r in 1:8, j in 1:30], x[r,j,"G\\/N"] <= x[r,j,"G"] + x[r,j,"N"]
    # no more than 2 consecutive G[N], G[N] 
    # TODO(odow): is this right? 
    [r in 1:8, g in 3:30], sum(x[r,i,"G\\/N"] for i in g-2:g) <= 2
    # no G can follow N 
    [r in 1:8, g in 2:30], x[r,g,"N"] <= 100 * (1 - x[r,g-1,"G"])
    [r in 1:8, g in 2:30], x[r,g,"G/\\N"] <= 100 * (1 - x[r,g-1,"G"])
end)
sol = [init_sol[r, g] ? "" : "X" for r in 1:8, g in 1:30]
for i in 1:8, j in 1:30
    if init_sol[i, j] == 0
        fix.(x[i, j, ["G", "N"]], 0; force = true)
    end
end
for (i, j) in init_sol_1
    fix(x[i, j, "G"], 0; force = true)
    sol[i, j] = "xg"
end
for (i, j) in init_sol_2
    fix(x[i, j, "N"], 0; force = true)
    sol[i, j] = "xn"
end
optimize!(model)
for i in 1:8, j in 1:30, k in ["G", "N"]
    if round(Int, value(x[i, j, k])) == 1
        sol[i, j] *= k
    end
end
for i in 1:8
    for j in 1:30
        print(lpad(sol[i, j], 4))
    end
    println()
end

rocco_sprmnt21 · May 28, 2023, 9:17pm

this should prevent sequences like:
G-G-G, G-N-G, N-GN-G, N-N-N, etc…

this should prevent sequences like:
N-G, N-GN, GN-G, GN-GN (so some constraints must be added)
But, I don’t understand why, it doesn’t seem to work for the set constraints either

odow · May 28, 2023, 9:30pm

I don’t really understand the details. But perhaps just:

    # no G can follow N 
    [r in 1:8, g in 2:30], x[r,g,"G"] <= 1 - x[r,g-1,"N"]

rocco_sprmnt21 · May 28, 2023, 9:38pm

tried to use this trick2

odow · May 28, 2023, 9:57pm

Because we know that all the variables are binary, the “big-M” can be 1. So remove the 100 * .

That gets your constraint to x[r,g,"N"] <= 1 - x[r,g-1,"G"]. Then swap sides: x[r,g-1,"G"] <= 1 - x[r,g,"N"].

This says that if an N shift in period g, then the G shift in period g-1 must be 0. I think that’s the wrong way around.

My constraint x[r,g,"G"] <= 1 - x[r,g-1,"N"] says that if N in period g-1, then the G in period g must be 0. So no G can follow N.

rocco_sprmnt21 · May 29, 2023, 7:14am

The fact that you can use “talking” indexes is very nice. It occurred to me but I didn’t think it was possible.
I would like to have someone confirm that the following extension of the translation of a multiple OR in linear constraints is correct.
Same for extending the constraints that “G” can’t come after “N”


    #  x["G\/N"] = x["G"] \/ x["N"]\/x["G/\\N"]
    [r in 1:8, j in 1:30], x[r,j,"G"] <= x[r,j,"G\\/N"]
    [r in 1:8, j in 1:30], x[r,j,"N"] <= x[r,j,"G\\/N"]
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] <= x[r,j,"G\\/N"]
    [r in 1:8, j in 1:30], x[r,j,"G\\/N"] <= x[r,j,"G"] + x[r,j,"N"]+ x[r,j,"G/\\N"]
    # no more than 2 elements consecutive from {G, N, GN} 
    [r in 1:8, g in 3:30], sum(x[r,i,"G\\/N"] for i in g-2:g) <= 2

    # no G can follow N 
    [r in 1:8, g in 2:30], x[r,g,"G"] <= 1 - x[r,g-1,"N"]
    [r in 1:8, g in 2:30], x[r,g,"G"] <= 1 - x[r,g-1,"G/\\N"]
    [r in 1:8, g in 2:30], x[r,g,"G/\\N"] <= 1 - x[r,g-1,"N"]
    [r in 1:8, g in 2:30], x[r,g,"G/\\N"] <= 1 - x[r,g-1,"G/\\N"]

odow · May 29, 2023, 9:32am

I would like to have someone confirm that the following extension of the translation of a multiple OR in linear constraints is correct.

Correct. It helps to walk through the logic, treating each variables 0 = false and 1+ = true

[r in 1:8, j in 1:30], x[r,j,"G"] <= x[r,j,"G\\/N"]

If LHS is 0, then RHS can be 0 or 1. But if LHS is 1, then RHS is 1. So this is G => G\/N

[r In 1:8, j in 1:30], x[r,j,"G\\/N"] <= x[r,j,"G"] + x[r,j,"N"] + x[r,j,"G/\\N"]

If LHS is 0, then RHS can be 0, 1, 2, or 3. If LHS is 1, then at least one of the RHS must be 1. Which is what you want.

rocco_sprmnt21 · May 29, 2023, 8:55pm

Let’s say that I get as input a table with constraints and boundary conditions (that is, already pre-compiled cells).
On these I apply the init_model() function to take into account all these conditions:

julia> m_
8×30 Matrix{Any}:
 "xx"  "xx"  "xx"  "xx"  "."   "."   "."   "."   "xx"  "."   "."   "."   "."   "."   …  "."   "."   "."   "."   "."   "."   "."   "."   "."   "."   "."   "G"   "N"      
 "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "."   "."   "."   "."   "."   "."   "."      "xx"  "xx"  "."   "."   "."   "."   "."   "."   "."   "."   "."   "xx"  "xx"     
 "G"   "xx"  "G"   "xg"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "."   "xx"  "xx"     "xg"  "."   "."   "xx"  "xx"  "xx"  "xx"  "."   "."   "."   "xx"  "xx"  "xx"     
 "xg"  "GN"  "N"   "xg"  "."   "xx"  "xx"  "."   "."   "xx"  "xx"  "."   "."   "."      "xx"  "."   "."   "."   "."   "."   "xx"  "xx"  "."   "."   "."   "N"   "."      
 "N"   "xx"  "xx"  "xx"  "."   "."   "xx"  "."   "xg"  "xg"  "."   "."   "."   "."      "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"     
 "xx"  "xx"  "xx"  "xx"  "xn"  "xg"  "xn"  "."   "xx"  "."   "."   "xx"  "."   "xn"  …  "."   "xx"  "."   "xn"  "."   "xx"  "xx"  "xx"  "xn"  "."   "xn"  "."   "xnG"
 "xx"  "xx"  "xx"  "GN"  "."   "."   "."   "."   "."   "."   "."   "."   "."   "."      "."   "."   "."   "."   "."   "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"  "xx"     
 "xx"  "xx"  "xx"  "xx"  "xx"  "."   "."   "."   "."   "."   "."   "."   "."   "."      "xx"  "."   "."   "."   "."   "xx"  "xx"  "xx"  "."   "."   "xx"  "xx"  "xx"  


function init_model(m)
    fixm(i,j, k,b)= begin fix(x[i, j, k], b; force = true); m[i,j]=replace(m[i,j], r"[N|G]"=>"") end
    for i in 1:8, j in 1:lastday
        m[i, j] == "xx"  ? fixm.(i, j, ["G", "N"], 0)       :
        m[i, j] == "xg"  ? fixm(i, j, "G", 0)               :
        m[i, j] == "xn"  ? fixm(i, j, "N", 0)               :
        m[i,j]∈["G","N"] ? fixm(i, j, m[i,j], 1)            :    
        m[i, j] == "xnG" ? fixm.(i, j, ["N","G"], [0,1])    :
        m[i, j] == "xgN" ? fixm.(i, j, ["G","N"], [0,1])    :
        m[i, j] == "GN"  ? fixm(i, j, "G/\\N", 1)           :
        nothing
    end
end

If I define the objective function as follows, is it guaranteed that the obtained solution is an absolute minimum?
It’s unique?

@objective(model, Min, sum(x[:,:,"G/\\N"]))

Summary

using XLSX, Dates

#path="xxxxxyyyyzzz.xlsx"
function getdata(path,sh,rngData,rngTot)
    xf = XLSX.readxlsx(path)
    data=xf[sh][rngData]
    tot=xf[sh][rngTot]
    m=replace(data, missing=>".")
    M=Dict("G" => tot[:,1], "N" => tot[:,2])
    m,M
end

function weekends(m,y)
    d=Date(Dates.Month(m), Dates.Year(y))
    day.(filter(d->dayofweek(d)==6||dayofweek(d)==7, d:lastdayofmonth(d)))
end

mese=6
lastday=day(lastdayofmonth(Date(Dates.Month(mese), Dates.Year(2023))))
H=(1,"B3:AE10","AI3:AJ10") #####  foglio excel,celle indispon. e tot per operatore
m, M = getdata(path,H[1],H[2],H[3])

m_=copy(m)
#----------------
using JuMP
import HiGHS

model = Model(HiGHS.Optimizer)
# set_silent(model)
@variable(model, x[i in 1:8, j in 1:30, k in ["G", "N", "G/\\N", "G\\/N"]], Bin)
@objective(model, Min, sum(x[:,:,"G/\\N"]))
@constraints(model, begin
    # total of G and N for row
    [k in ["G", "N"], r in 1:8], sum(x[r,:,k]) == M[k][r]
    # each day must have exactly one G (k==1) and one N (k==2)
    [k in ["G", "N"], g in 1:30], sum(x[:,g,k]) == 1
    # max of G and N for row for weekends
    #[k in ["G", "N", "G/\\N"], r in 1:8], sum(x[r,j,k] for j in wse) <= 3
    [ r in 1:8], sum(x[r,j,"G\\/N"] for j in weekends(mese,2023)) <= 3
    # max one GN for row
    [r in 1:8], sum(x[r,g,"G/\\N"] for g in 1:30) <= 1
    #  x["G/\N"] = x["G"] /\ x["N"]
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] <= x[r,j,"G"]
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] <= x[r,j,"N"]           
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] >= x[r,j,"G"] + x[r,j,"N"] - 1

    #  x["G\/N"] = x["G"] \/ x["N"]\/x["G/\\N"]
    [r in 1:8, j in 1:30], x[r,j,"G"] <= x[r,j,"G\\/N"]
    [r in 1:8, j in 1:30], x[r,j,"N"] <= x[r,j,"G\\/N"]
    [r in 1:8, j in 1:30], x[r,j,"G/\\N"] <= x[r,j,"G\\/N"]
    [r in 1:8, j in 1:30], x[r,j,"G\\/N"] <= x[r,j,"G"] + x[r,j,"N"]+ x[r,j,"G/\\N"]
    # no more than 2 elements consecutive from {G, N, GN} 
    [r in 1:8, g in 3:30], sum(x[r,i,"G\\/N"] for i in g-2:g) <= 2
    # no G can follow N 
    [r in 1:8, g in 2:30], x[r,g,"G"] <= 1 - x[r,g-1,"N"]
    [r in 1:8, g in 2:30], x[r,g,"G"] <= 1 - x[r,g-1,"G/\\N"]
    [r in 1:8, g in 2:30], x[r,g,"G/\\N"] <= 1 - x[r,g-1,"N"]
    [r in 1:8, g in 2:30], x[r,g,"G/\\N"] <= 1 - x[r,g-1,"G/\\N"]
end)

function init_model(m)
    fixm(i,j, k,b)= begin fix(x[i, j, k], b; force = true); m[i,j]=replace(m[i,j], r"[N|G]"=>"") end
    for i in 1:8, j in 1:lastday
        m[i, j] == "xx"  ? fixm.(i, j, ["G", "N"], 0)       :
        m[i, j] == "xg"  ? fixm(i, j, "G", 0)               :
        m[i, j] == "xn"  ? fixm(i, j, "N", 0)               :
        m[i,j]∈["G","N"] ? fixm(i, j, m[i,j], 1)            :    
        m[i, j] == "xnG" ? fixm.(i, j, ["N","G"], [0,1])    :
        m[i, j] == "xgN" ? fixm.(i, j, ["G","N"], [0,1])    :
        m[i, j] == "GN"  ? fixm(i, j, "G/\\N", 1)           :
        nothing
    end
end

init_model(m)

optimize!(model)

for i in 1:8, j in 1:30, k in ["G", "N"]
    if round(Int, value(x[i, j, k])) == 1
        m[i, j] *= k
    end
end


XLSX.openxlsx(path, mode="rw") do xf
    sheet = XLSX.addsheet!(xf, "alt-19")
    XLSX.writetable!(sheet, DataFrame(m,:auto), anchor_cell=XLSX.CellRef("B2"))  
end


for i in 1:8
    for j in 1:30
        print(lpad(m[i, j], 4))
    end
    println()
end

odow · May 29, 2023, 10:22pm

If I define the objective function as follows, is it guaranteed that the obtained solution is an absolute minimum?

Yes

It’s unique?

No. There might be other solutions with the same objective value.