I would like to know if there is a synthetic form to express the constraint Max(Si)-Min(Si) <=s, without having to explicitly write all the possible differences (Si-Sj) and impose the constraint on these
S1:: x[1,1,N]+x[1,2,N]+ … +x[1,30,N]+ x[1,1,G]+x[1,2,G]+ … +x[1,30,G]
S2:: x[2,1,N]+x[2,2,N]+ … +x[2,30,N]+ x[2,1,G]+x[2,2,G]+ … +x[2,30,G]
…
S5:: x[5,1,N]+x[5,2,N]+ … +x[5,30,N]+ x[5,1,G]+x[5,2,G]+ … +x[5,30,G]
Max(Si)-Min(Si) <=s
Perhaps, have a collection of constraints:
(adding a couple of variables m and M)
M - Si >= 0 for all i
Si - m >= 0 for all i
M >= 0
m >= 0
M - m <= s
I’m not sure I understand the implementation detail.
My variables are all in {0,1}.
How do I express M and m?
Make m and M integers. (The sums S_i are integers potentially).
This is the current model.
nr= 5 (but it is an aprameter that sometimes changes)
ng= 30 (in general these are the days of the current month).
Can I mix binary variables and integer variables in the same model?
@variable(model, x[i in 1:nr, j in 1:ng, k in ["G", "N", "G/\\N", "G\\/N"]], Bin)
#ottimizza cercando il minimo GN per tutti tranne (xxxx e yyyy)per cui cerca il massimo numero di GN
@objective(model, Min, [sum(x[[1,2,5],1:ng,"G/\\N"]), - sum(x[[3,4],1:ng,"G/\\N"])])
You can in the model. The specific solver needs to handle it. But perhaps an MWE would help (i.e. something paste-able). And also odow might chime in.
I’ll try. I hope I don’t make a mess.
I also TRY TO attach an excel file to support the example.
For each day (matrix column) resources must be allocated to cover service N and service G. They can be two distinct ops or the same op.
There are two ops that prefer to have as many GNs as possible (i.e. double services assigned to the same op on the day) while for the others the objective is to minimize double services (GNs) on the same day.
The input matrix (built on the basis of an Excel table) provides for partial or total unavailability that the various operators express for the current month and there are also preliminary assignments of services for some operators on some days.
In the attached excel table there is a legend that explains the meaning of the symbols used.
The constraints used tend to construct a “balanced” allocation solution.
The purpose of adding the constraint in question is to improve this balance.
I would like to give it a try and look for a solution that shortens the distance between the max(si) and min(Si). In the case of the example max(Si)=7+6 for op2 and min(Si)=5+5 for op5.
If it exists, I would like to look for a solution that, all other constraints being equal(*), limits this distance to 2.
It is required that no G can follow N for the same operator. and others (which I can illustrate if necessary).
Among the constraints I added the “tricks” relating to the logical conditions OR and AND
(*) for example, it is required not to have too many consecutive allocations or too many allocations on weekends and/or Sundays.
using XLSX, Dates
using CSV, Tables, DataFrames
path="febbraio 2024 MWE.xlsx"
function getdata(path,sh,rngData,rngTot)
xf = XLSX.readxlsx(path)
#sn=XLSX.sheetnames(xf)
data=xf[sh][rngData]
mM=xf[sh][rngTot]
m=replace(data, missing=>".")
RmM=Dict("min" => mM[:,1], "max" => mM[:,2])
m,RmM
end
function weekends(m,y)
#we = [6, 7, 13, 14, 20, 21, 27, 28, 34, 35]
d=Date(Dates.Month(m), Dates.Year(y))
#lastday=day(lastdayofmonth(d))
#startday=(firstdayofmonth(d) - firstdayofweek(d)).value
# startday=dayofweek(d)-1
# filter(d-> 0 < d < lastday+1, we .- startday)
day.(filter(d->dayofweek(d)==6||dayofweek(d)==7, d:lastdayofmonth(d)))
#day.(filter(d->Dates.issaturday(d)||Dates.issunday(d), d:lastdayofmonth(d)))
end
function sundays(m,y)
# d=Date(Dates.Month(m), Dates.Year(y))
# day.(filter(d->dayofweek(d)==7, d:lastdayofmonth(d)))
weekends(m,y)[2:2:end]
end
function holidays(m,y)
d=Date(Dates.Month(m), Dates.Year(y))
day.(filter(d->dayofweek(d)==7, d:lastdayofmonth(d)-Day(1)))
end
nr=5
mese=2
y=2024
ng=daysinmonth(Date(y,mese))
col='A'*('A'+ng-26)*"$(nr+2)" # last cell
H=(1,"B3:$col","AI3:AJ$(nr+2)") ##### foglio excel,celle indispon. e min max for op
m, RmM = getdata(path,H[1],H[2],H[3])
#----------------
using JuMP
import HiGHS
import MultiObjectiveAlgorithms as MOA
model = Model()
set_silent(model)
@variable(model, x[i in 1:nr, j in 1:ng, k in ["G", "N", "G/\\N", "G\\/N"]], Bin)
@objective(model, Min, [sum(x[[1,2,5],1:ng,"G/\\N"]), - sum(x[[3,4],1:ng,"G/\\N"])])
@constraints(model, begin
# min max of G and N for row
[k in ["G","N"], r in 1:nr], sum(x[r,:,k]) >= RmM["min"][r]
[k in ["G","N"], r in 1:nr], sum(x[r,:,k]) <= RmM["max"][r]
# sum(x[:,:,"G"]) == ng; # sum(x[:,:,"N"]) == ng questa é implicatadalla seguente # each day must have exactly one G and one N
daycom[k in ["G", "N"], g in 1:ng], sum(x[:,g,k]) == 1
# [r in 1:nr] max(M["G"][r]+M["N"][r]) - min(M["G"][r]+M["N"][r]) <=2 ????
# max num of G and N for row for sundays
sndlim[r in 1:nr], sum(x[r,j,"G"]+x[r,j,"N"]+2*x[r,j,"G/\\N"] for j in sundays(mese,y)) <= 5
#min and max num of G and N for row for weekends eccetto per op4
minwklim[r in 1:nr], sum(x[r,j,"G"]+x[r,j,"N"] for j in weekends(mese,y)) >= 3
maxwklim[r in [1,2,3,5]], sum(x[r,j,"G"]+x[r,j,"N"] for j in weekends(mese,y)) <= 5
# min and max of N for row for weekends eccetto per op4
[ r in [1,2,3,5]], sum(x[r,j,"N"] for j in weekends(mese,y)) >= 1
[ r in [1,2,3,5]], sum(x[r,j,"N"] for j in weekends(mese,y)) <= 3
# max 2 GN for row (escluso op3 and op4)
gnlim[r in [1,2,5]], sum(x[r,g,"G/\\N"] for g in 1:ng) <= 2
# # x["G/\N"] = x["G"] /\ x["N"]
[r in 1:nr, j in 1:ng], x[r,j,"G/\\N"] <= x[r,j,"G"]
[r in 1:nr, j in 1:ng], x[r,j,"G/\\N"] <= x[r,j,"N"]
[r in 1:nr, j in 1:ng], 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:nr, j in 1:ng], x[r,j,"G"] <= x[r,j,"G\\/N"]
[r in 1:nr, j in 1:ng], x[r,j,"N"] <= x[r,j,"G\\/N"]
#[r in 1:nr, j in 1:ng], x[r,j,"G/\\N"] <= x[r,j,"G\\/N"]
[r in 1:nr, j in 1:ng], x[r,j,"G\\/N"] <= x[r,j,"G"] + x[r,j,"N"] # + x[r,j,"G/\\N"]
# no more than 3 elements consecutive from {G, N, GN}
[r in 1:nr, g in 4:ng], sum(x[r,i,"G\\/N"] for i in g-3:g-1) <= 2
#[r in 1:nr, g in 4:ng], sum(x[r,j,"G"]+x[r,j,"N"]+2*x[r,j,"G/\\N"] for j in g-3:g-1) <= 2
# no 2 elements GN consecutive
[r in 1:nr, g in 2:ng], x[r,g,"G/\\N"] <= 1 - x[r,g-1,"G/\\N"]
# no G can follow N
[r in 1:nr, g in 2:ng], x[r,g,"G"] <= 1 - x[r,g-1,"N"]
[r in 1:nr, g in 2:ng], x[r,g,"G"] <= 1 - x[r,g-1,"G/\\N"]
[r in 1:nr, g in 2:ng], x[r,g,"G/\\N"] <= 1 - x[r,g-1,"N"]
[r in 1:nr, g in 2:ng], x[r,g,"G/\\N"] <= 1 - x[r,g-1,"G/\\N"]
# no G in sunday if N in saturday (se sabato(N) => !G(domenica))
#[r in 1:nr, g in weekends(mese,y)[1:2:end]], x[r,g+1,"G"] <= 1 - x[r,g,"N"]
#[r in 1:nr, g in holidays(mese,y)], x[r,g+1,"G/\\N"] <= 1 - x[r,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:nr, j in 1:ng
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)
set_optimizer(model, () -> MOA.Optimizer(HiGHS.Optimizer))
set_attribute(model, MOA.Algorithm(), MOA.EpsilonConstraint())
optimize!(model)
solution_summary(model)
mres=[]
for res in 1:result_count(model)
m_=copy(m)
for i in 1:nr, j in 1:ng, k in ["G", "N"]
if round(Int, value(x[i, j, k];result=res)) == 1
m[i, j] *= k
end
end
push!(mres,m)
m=copy(m_)
end
function sinossi(m)
nGN=zeros(Int,nr,6)
for r in 1:nr
nGN[r,1]=count(contains("G"),m[r,weekends(mese,y)[1:2:end]])
nGN[r,2]=count(contains("N"),m[r,weekends(mese,y)[1:2:end]])
nGN[r,3]=count(contains("G"),m[r,weekends(mese,y)[2:2:end]])
nGN[r,4]=count(contains("N"),m[r,weekends(mese,y)[2:2:end]])
nGN[r,5]=count(contains("G"),m[r,:])
nGN[r,6]=count(contains("N"),m[r,:])
end
nGN
end
#----------------------
# using CSV, Tables, DataFrames
#CSV.write("turni_giugno_sol_jump.csv", Tables.table(m),delim=',')
#--------------------
for res in 1:result_count(model)
XLSX.openxlsx(path, mode="rw") do xf
sheet = XLSX.addsheet!(xf, "result_$(res)_xc2")
XLSX.writetable!(sheet, DataFrame(mres[res],:auto), anchor_cell=XLSX.CellRef("B2"))
XLSX.writetable!(sheet, DataFrame(sinossi(mres[res]),["GS","NS","GD","ND","totG","totN"]), anchor_cell=XLSX.CellRef("AH2"))
end
end
This example is not so minimal… in any case, I’ve made some example which might clarify my suggestion:
julia> model = Model(HiGHS.Optimizer);
julia> set_silent(model);
julia> @variable(model, x[1:10], Bin)
julia> using SparseArrays
julia> sp = sprand(Bool,20,10,0.5).*(-1).^rand(Bool, 20,10)
20×10 SparseMatrixCSC{Int64, Int64} with 105 stored entries:
⎡⠣⣲⣻⢞⡦⎤
⎢⠎⡤⣪⠐⠢⎥
⎢⠫⠦⣧⡑⡨⎥
⎢⣩⣧⢏⢦⣋⎥
⎣⡿⠟⠃⡪⣋⎦
julia> cvec = sp*x .+ rand(0:10,20)
20-element Vector{AffExpr}:
x[1] + x[5] + x[6] + x[8] + 3
-x[1] - x[3] - x[4] + x[5] - x[6] + x[7] + x[8] - x[9] + 4
x[2] - x[4] + x[6] + x[7] + x[9] + x[10] + 5
x[3] - x[4] + x[5] + x[6] - x[8] - x[9] + 2
-x[2] + x[6] + 6
-x[1] - x[5] + x[8] - x[9] + 7
-x[1] - x[3] + x[4] + x[6] + x[10] + 3
x[3] + x[5] - x[6] + 6
x[1] - x[2] + x[5] + x[7] - x[10] + 10
x[1] + x[3] + x[5] + x[8] + 8
x[2] - x[3] + x[4] - x[5] + x[6] - x[10] + 6
-x[5] + x[6] - x[7] + x[9] + 8
x[1] + x[2] - x[3] + x[5] - x[6] + x[9] + x[10] + 10
x[3] - x[5] + x[7] - x[9] + 3
x[2] - x[3] + x[4] + x[5] - x[7] + x[8] + 6
x[1] + x[2] - x[3] - x[4] + x[6] + x[8] + x[9] - x[10] + 3
-x[1] - x[2] + x[3] + x[4] - x[5] - x[8] + x[9] + x[10] + 9
x[1] - x[2] - x[3] + x[4] + x[5] - x[7] + x[9]
-x[1] - x[2] + x[3] - x[8]
-x[1] + x[7] - x[9] - x[10] + 3
julia> @variable(model, M);
julia> @variable(model, m);
julia> for i in eachindex(cvec)
@constraint(model, cvec[i] <= M)
@constraint(model, cvec[i] >= m)
end
julia> @objective(model, Min, M - m)
M - m
julia> optimize!(model)
julia> solution_summary(model)
* Solver : HiGHS
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"kHighsModelStatusOptimal"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : NO_SOLUTION
Objective value : 1.00000e+01
Objective bound : 1.00000e+01
Relative gap : 0.00000e+00
* Work counters
Solve time (sec) : 2.59709e-03
Simplex iterations : 5
Barrier iterations : -1
Node count : 1
julia> value.(x)
10-element Vector{Float64}:
0.0
0.0
1.0
0.0
1.0
1.0
0.0
0.0
0.0
1.0
julia> value(M)
10.0
julia> value(m)
0.0
So the M and the m bound the differnece in expression values, but need not be {0,1}, in fact, they can even be real. If there is a specific range between maximum and minimum, this can be encoded with a @constraint M - m <= s.
Hi @Dan and @rocco_sprmnt21, @Dan’s post above is that I would suggest.
I would perhaps even go as far as testing:
@variable(model, cvec[1:20])
@constraint(model, cvec .== sp*x .+ rand(0:10,20))
Sometimes a sparser formulation can be more efficient, even if it means adding extra variables.
If you find the time and inclination to take a look at my script (not minimal at all, I must admit) can you tell me whether they are useful and, if so, how I can use sparse matrix forms?
I would also like to know if there is a way in calculating the optimum subject to a list of constraints to have a suboptimal result by “deactivating” the last (or last k) constraints (assuming that the constraints are used according to the order in which they appear. But maybe this isn’t the case)?
Obviously I can do this operation by hand. but… it would be useful for me to have an automatic, even semi-automatic mode 
in other words, are “elastic constraints” possible?