Non-convergent "weighted" objective function

Specific Domains Optimization (Mathematical)
1

I need to find a way to allocate a pool of resources to perform a series of daily services over the course of a month.
There are obviously a long series of general (for all resources) and specific (for individual resources) constraints.
The objective function to be minimized tries to make the assignment “balanced.”
Using HIGHS, I can’t use quadratic functions, but I approximate the result by adding two auxiliary variables (lower_bound=m and upper_bound=M) and imposing Min(M-m).
The script converges quickly and provides a result.

...
model = Model(HiGHS.Optimizer)
#set_silent(model)
@variable(model, x[1:nr, 1:ng, turni ], Bin)
@variable(model, z[1:nr, 1:ng], Bin)
@constraints(model, begin

    # each day must have exactly one Ti or one Fi
    fer1[k in union(t_feriali,t_card), g in g_feriali], sum(x[:,g,k]) == 1
    fer2[k in t_festivi, g in g_feriali], sum(x[:,g,k]) == 0

    fest1[k in t_festivi, g in g_festivi], sum(x[:,g,k]) == 1
    fest2[k in union(t_feriali, t_card), g in g_festivi], sum(x[:,g,k]) == 0

    sat[k in t_Cs, g in saturdays(mese,y)], sum(x[:,g,k]) == 1
    no_sat[k in ["CUMs","CLMs"], g in setdiff(1:ng,saturdays(mese,y))], sum(x[:,g,k]) == 0
...

    # no more than 5 turni consecutive  
    [r in 1:nr, g in 6:ng], sum(sum(x[r,i,k] for k in turni ) for i in g-5:g) <= 5


    # nothing can follow N or NF
     [r in 1:nr, g in 2:ng, k in turni], x[r,g,k] <= 1 - x[r,g-1,"NF"]
     [r in 1:nr, g in 2:ng, k in turni], x[r,g,k] <= 1 - x[r,g-1,"N"]

    #distance between N and N or NF >=4
     [r in 1:nr, g in 5:ng], sum(x[r,i,"N"]+x[r,i,"NF"] for i in g-4:g) <= 1
     
end)


peso = Dict{Int, Float64}()
for g in g_festivi
        peso[g] = 1.5
end
for g in g_feriali
        peso[g] = 1.0
end


# Carico di lavoro
@expression(model, workload[r=no_E], sum(x[r, g,k]*peso[g] for g in 1:ng, k in turni))

# Variabili per max/min carico
@variable(model, M)   # massimo carico
@variable(model, m)   # minimo carico

@constraint(model, [r in no_E], workload[r] <= M)
@constraint(model, [r in no_E], workload[r] >= m)

# Obiettivo: minimizzare squilibrio massimo
@objective(model, Min, M - m)
...

To obtain a more “fair” result, I tried solving the same problem by assigning different weights to some days of the month.
In this form, I obtain an “infinite” sequence of the following lines:

funzione obiettivo non converge
funzione obiettivo non converge1107×407 14.9 KB

How can I read this output?
Can you give me some general suggestions for overcoming this problem?
If needed, I can post the complete script.

2

An MIP formulation may become harder to solve if you modify the coefficient of the objective expression. In your case it seems that the solver is no longer making progress. And the Gap is so huge… I think you need to double check your formulation. For a problem in practice, how can you get a 60% gap…

3

I have no idea what the GAP parameter means.
I tried changing the weights to only integer values
1 and 2 instead of 1 and 1.5, and it converged quickly.
PS
Is there a way to stop the calculation afterwards without aborting the terminal?
I mean using some appropriate parameter in the optimize!() function call.

Or similarly at a higher level.

PPS
It seems to not handle non-integer numbers.

I only have two sets of parameters.
If I set s1=1 and s2=1.5, it doesn’t converge.
If I set s1=2 and s2=3, it converges.

My implication is that the two formulations are mathematically equivalent.

4

ctrl + c is user’s interrupt, which is safe—you can re-call optimize! later.

set some attributes before you optimize!, e.g. TimeLimit, NodeLimit, …

You need to have some basic knowledge about global optimization, e.g. via reading books, articles, or tutorials…

At least for Gurobi, I think its API only accepts data in Float64 form, so they should have no difference. (I don’t know how HiGHS does)

Finally, I think you may consider revise your first post and give the runnable code so others can understand your problem better.

5

Here’s the script; unfortunately, I have no way to simplify it.
Lines 180-190 contain the ones that change the algorithm’s behavior.

Summary 
using XLSX, Dates , CSV, Tables, DataFrames, JuMP,  HiGHS

function sol_jump(mese,y)

function getdata(path,sh,rngData)
    xf = XLSX.readxlsx(path)
    #sn=XLSX.sheetnames(xf)
    data=xf[sh][rngData]
    replace(data, missing=>".")
end

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

function mondays(m,y)
    d=Date(Dates.Month(m), Dates.Year(y))
    day.(filter(d->dayofweek(d)==1, d:lastdayofmonth(d)))
end

 function MNDS(m,y)
     d=Date(Dates.Month(m), Dates.Year(y))
     starts=max(d,firstdayofweek(d))
     mondays=filter(d->dayofweek(d)==1, d:lastdayofmonth(d))
     ends=tonext.(d->dayofweek(d) == Dates.Friday, mondays)
     ends[end]=min(ends[end], lastdayofmonth(d))
     (:).(day.(mondays),day.(ends))
 end
 wks=MNDS(mese,y)
#controllare se funziona per i diversi mesi.
# se nel mese c'è prima un venerdì di un lunedì?
# se lunedi+ 4 supera 31?
# ecc

ops=["BAT",
"CAVA",
"CEL",  
"DE_V",  
"EMO",
"GAT",
"LAN",
"LUO",  
"MAL",
"MAN",
"PAG",
"RIZ",
"SOL"
]


BAT = findfirst(==( "BAT"), ops,) # 1
CEL = findfirst(==( "CEL"), ops,) # 3
DE_V = findfirst(==( "DE_V"), ops,) # 4
EMO = findfirst(==( "EMO"), ops,) # 5
GAT = findfirst(==( "GAT"), ops,) # 6
LUO = findfirst(==( "LUO"), ops,) # 8
MAN = findfirst(==( "MAN"), ops,) # 10
PAG = findfirst(==( "PAG"), ops,) # 11

nr = length(ops)
ops_lim=setdiff(1:nr,[1,3,4,5,8,11])
no_E=setdiff(1:nr,EMO)
ops_no_CUMs_CLMs=setdiff(1:nr,[4,8,11])  # tutti meno  DE_V, LUO, PAG

ng = daysinmonth(Date(y,mese))

col='A'*('A'+ng-26)*"$(nr+2)" # colonna AE10 o AGF0
H = (1,"B3:$col") #####  foglio excel,celle indispon. e ranges per persona
res= getdata(path,H[1],H[2])

g_festivi = weekends(mese,y)
g_feriali = setdiff(1:ng,g_festivi)

turni = ["M1","M2", "P", "N", "GF","NF","CUM1","CUM2","CLM","CUMs","CLMs"]
t_feriali = turni[1:4]
t_festivi = turni[5:6]
t_card = turni[7:9]
t_Cs = turni[10:11]
t_mattina = ["M1","M2","CUM1","CUM2","CLM","GF","CUMs","CLMs"]

model = Model(HiGHS.Optimizer)
#set_silent(model)
@variable(model, x[1:nr, 1:ng, turni ], Bin)
@variable(model, z[1:nr, 1:ng], Bin)
@constraints(model, begin

    # each day must have exactly one t: (t in turni)
    fer1[k in union(t_feriali,t_card), g in g_feriali], sum(x[:,g,k]) == 1
    fer2[k in t_festivi, g in g_feriali], sum(x[:,g,k]) == 0

    fest1[k in t_festivi, g in g_festivi], sum(x[:,g,k]) == 1
    fest2[k in union(t_feriali, t_card), g in g_festivi], sum(x[:,g,k]) == 0

    sat[k in t_Cs, g in saturdays(mese,y)], sum(x[:,g,k]) == 1
    no_sat[k in ["CUMs","CLMs"], g in setdiff(1:ng,saturdays(mese,y))], sum(x[:,g,k]) == 0

    # range of N and NF for row for month
    [r in ops_lim], sum(x[r,j,"N"]+x[r,j,"NF"] for j in 1:ng) <= 4  #4
    [r in ops_lim], sum(x[r,j,"N"]+x[r,j,"NF"] for j in 1:ng) >= 3  #3
    
    [r in ops_lim], sum(x[r,j,"N"] for j in 1:ng) <= 3
    [r in ops_lim], sum(x[r,j,"N"] for j in 1:ng) >= 1
   
    # max of  NF for row for month
    [ r in 1:nr], sum(x[r,j,"NF"] for j in g_festivi) <= 2  # 2
    
    # NF==2 --> N<=1
    [r in ops_lim], sum(x[r,:,"NF"]) == sum(i * z[r, i] for i in 1:ng)
    [r in ops_lim], z[r,2] --> {sum(x[r, :, "N"]) <= 1}
    [r in ops_lim], sum(z[r,:]) <= 1


    [r in 1:nr, w in wks], x[r,w[1],"CUM1"] --> {sum(x[r,Base.rest(w),"CUM1"])>=length(Base.rest(w))}
    [r in 1:nr, w in wks], x[r,w[1],"CUM2"] --> {sum(x[r,Base.rest(w),"CUM2"])>=length(Base.rest(w))}
    [r in 1:nr, w in wks], x[r,w[1],"CLM" ] --> {sum(x[r,Base.rest(w),"CLM" ])>=length(Base.rest(w))}

    # max turni nello stesso giorno per operatore
    # non notte insieme a mattina o pomeriggio, ecc.
    max_turni_fer_N_P[op in 1:nr, g in g_feriali], sum(x[op,g,k]  for k in ["N","P"]) <= 1
    max_turni_fer_N_M[op in 1:nr, g in g_feriali], sum(x[op,g,k]  for k in ["N";t_mattina]) <= 1
    max_turni_fer_P_M[op in 1:nr, g in g_feriali], sum(x[op,g,k]  for k in ["P";t_mattina]) <= 1
    max_turni_fer[op in 1:nr, g in g_feriali], sum(x[op,g,k]  for k in t_feriali) <= 2
    max_turni_CU[op in 1:nr, g in g_feriali], sum(x[op,g,k]  for k in t_card) <= 1
    max_turni_fest[op in 1:nr, g in g_festivi], sum(x[op,g,k]  for k in [t_festivi;t_Cs]) <= 1

    # max num of GF and NF for row for weekends
    wk_lim[r in 1:nr], sum(x[r,j,"GF"]+x[r,j,"NF"] for j in weekends(mese,y)) <= 2  #2
  
#     # range for month # sostituiti dalla funzione obiettivo
#     [r in no_E], sum(sum(x[r,j,k] for k in turni ) for j in 1:ng) <= 16  #16
#     [r in no_E], sum(sum(x[r,j,k] for k in turni ) for j in 1:ng) >= 14  #14
    
    # vincoli sui turni specifici per operatore
    # es. BAT 3 turni di notte a mese
    # PAG, LUO, CEL MAN NON COPRONI I TURNI UTIC DI MATTINA
    # BAT SEMPRE IN UTIC LA MATTINA TRANNE SE FA LA NOTTE

    sum(x[BAT,g,"N"] + x[BAT,g,"NF"] for g in 1:ng) == 3     # 3 notti
    sum(x[BAT,g,"P"]  for g in 1:ng)== 0                             # no pomeriggi
    sum(x[DE_V,g,k] for k in ["N", "NF"],  g in 1:ng) == 0            # no notti
    sum(x[LUO,g, k] for k in ["N", "NF","M1","M2"],  g in 1:ng) == 0      # no notti, no UTIC mattina
    4 <= sum(x[LUO,g,"P"] for g in 1:ng) <= 5                             # 4 o 5 pomeriggi
    sum(x[CEL,g,k] for k in ["GF", "NF","M1","M2"], g in 1:ng) == 0      # no festivi, no UTIC mattina
    3 <= sum(x[CEL,g,"P"] for g in 1:ng) <= 4                            # 3 o 4 pomeriggi
    sum(x[CEL,g,"N"] for g in 1:ng) ==2                                  # 2 notti
    sum(x[CEL,g,k] for k in ["M1","M2"], g in 1:ng) == 0                 # no UTIC mattina
    sum(x[MAN,g,k] for k in ["M1","M2"], g in 1:ng) == 0                # no UTIC mattina
    4 <= sum(x[PAG,g,"P"] for g in 1:ng) <= 5                            # 4 o 5 pomeriggi
    sum(x[PAG,g,k] for k in ["N","GF","NF","M1","M2"] , g in 1:ng) == 0  # no notti, no festivi, no UTIC mattina
    sum(x[GAT,g,t]  for g in 1:ng, t in t_card) == 0                       # no cardio urgenza

    # prevalentemente DE_V, LUO, PAG  dovrebbero fare CU di sabato

    sum(x[r,g,"CUMs"]+x[r,g,"CLMs"] for g in saturdays(mese,y),r in ops_no_CUMs_CLMs) <= 1

    # no more than 5 turni consecutive  
    [r in 1:nr, g in 6:ng], sum(sum(x[r,i,k] for k in turni ) for i in g-5:g) <= 5

    # nothing can follow N or NF
     [r in 1:nr, g in 2:ng, k in turni], x[r,g,k] <= 1 - x[r,g-1,"NF"]
     [r in 1:nr, g in 2:ng, k in turni], x[r,g,k] <= 1 - x[r,g-1,"N"]

    #distance between N and N or NF >=4
     [r in 1:nr, g in 5:ng], sum(x[r,i,"N"]+x[r,i,"NF"] for i in g-4:g) <= 1
     
end)

peso = Dict{String, Int}()  #  THIS IS CONVERGENT
for g in turni
     if g in [t_festivi; "N"]
        peso[g] = 2
     else
        peso[g] = 1
     end
end



# peso = Dict{String, Float64}() THIS DOES NOT CONVERGE
# for g in turni
#      if g in [t_festivi; "N"]
#         peso[g] = 1
#      else
#         peso[g] = 0.5
#      end
# end



# Carico di lavoro
@expression(model, workload[r=no_E], sum(x[r, g,k]*peso[k] for g in 1:ng, k in turni))

# Variabili per max/min carico
@variable(model, M)   # massimo carico
@variable(model, m)   # minimo carico

@constraint(model, [r in no_E], workload[r] <= M)
@constraint(model, [r in no_E], workload[r] >= m)

# Obiettivo: minimizzare squilibrio massimo
@objective(model, Min, M - m)

function init_model(res)
    fixm(i,j, k,b)= begin fix(x[i, j, k], b; force = true); res[i,j]=replace(res[i,j], r"[M|P|N|E]"=>".") end
    for i in 1:nr, j in 1:ng
        res[i, j] == "X"  ? fixm.(i, j, turni, 0)                          :
        res[i, j] == "M"  ? fixm.(i, j, t_mattina, 0)                      :
        # chiarire se res[i, j] == "P" esclude anche "GF"                            #
        res[i, j] == "P"  ? fixm.(i, j, ["P"], 0)                          :
        res[i, j] == "N"  ? fixm.(i, j, ["N","NF"], 0)                     :  
        res[i, j] == "MP" ? fixm.(i, j, ["P";t_mattina], 0)                :
        res[i, j] == "PN" ? fixm.(i, j, ["P","N","NF"], 0)                 :
        res[i, j] == "MN" ? fixm.(i, j, ["N";t_mattina], 0)                :
        res[i, j] == "E" ? fixm.(i, j, turni, [0,0,0,1,0,0,0,0,0,0,0])     :
        nothing
    end
    for i in 1:nr, j in 1:ng
        if (res[i,j]==".."||res[i,j]=="...") 
            res[i,j]="."
        end
    end
end

init_model(res)
optimize!(model)
assert_is_solved_and_feasible(model)
x,nr,ng,turni,ops,res
end  # soljump



function report_to_excel(path, x,nr,ng,turni,ops,res,shn="1")

for i in 1:nr, j in 1:ng, k in turni
    if value(x[i, j, k]) >= 0.5
        if res[i,j]=="."    
            res[i, j] = k
        else
            res[i, j] *= "+"*k
        end
    end
end

vv=["" for i in 1:ng, j in 1:length(turni)]
resp=DataFrame(vv,turni)
for  j in 1:ng, k in turni, i in 1:nr
    if value(x[i, j, k]) >= 0.5
        resp[j, k] *= ops[i]
    end
end

giormi_it=Dict(
  "Wed" => "Mer",
  "Tue" => "Mar",
  "Thu" => "Gio",
  "Sun" => "Dom",
  "Mon" => "Lun",
  "Fri" => "Ven",
  "Sat" => "Sab")

git=[giormi_it[d] for d in Dates.format.(Date.(y, mese, 1:ng),"e")]
gn = Dates.format.(Date.(y, mese, 1:ng),"dd")
insertcols!(resp,1, :giorno => gn .* "-" .* git)


# output mode vista operatore
XLSX.openxlsx(path, mode="rw") do xf
    sheet = XLSX.addsheet!(xf, "sol_jump_m_"*shn*"p")
    XLSX.writetable!(sheet, resp, anchor_cell=XLSX.CellRef("A1"))  
end

# output mode vista servizio
XLSX.openxlsx(path, mode="rw") do xf
    sheet = XLSX.addsheet!(xf, "sol_jump_"*shn)
    # dnames = Dates.format.(Date.(y, mese, 1:ng),"dd") .* "-" .* (giormi_it[d] for d in Dates.format.(Date.(y, mese, 1:ng),"e"))
    dnames =  gn .* "-" .* git
    sol=DataFrame(res,dnames)
    insertcols!(sol,1, "giorno\\operatore" => ops)
    XLSX.writetable!(sheet, sol, anchor_cell=XLSX.CellRef("A1"))  
end


# output mode vista sintesi/sommario
XLSX.openxlsx(path, mode="rw") do xf
    sheet = XLSX.addsheet!(xf, "summary_jump_"*shn)
    vx=round.(Int,value.(x))
    df=DataFrame([ t=>[sum(vx[j, :, t]) for j in 1:nr] for t in turni])
    transform!(df,Cols(:) =>ByRow((x...)->sum(x)) => "tot")
    transform!(df,["N","NF"]=>ByRow((x...)->sum(x)) => "tot_N")
    insertcols!(df,1, :operatore => ops)
    XLSX.writetable!(sheet, df, anchor_cell=XLSX.CellRef("A1"))  
end

end


y, mese = 2025,10

path="C:\\Users\\sprmn\\OneDrive\\Desktop\\Turni San Giovanni\\test_ott_2025.xlsx"

x,nr,ng,turni,ops,res=sol_jump(mese,y)

report_to_excel(path, x,nr,ng,turni,ops,res,"peso_turni_1x2")

Here is the Excel table that the script uses to load the boundary conditions

test_ott_2025.xlsx

In the same Excel folder the script loads the processing results, in 3 separate sheets

6

Hi @rocco_sprmnt21,

There are three important columns to look at (you can ignore the rest):

  1. BestSol: this is the objective value of the best integer solution HiGHS has found so far. Your current best objective is 0.5.
  2. BestBound: this is a lower bound on the optimal objective value. So far, HiGHS has proven that there cannot be a solution with a cost lower than 0.197.
  3. Gap: this is the relative difference between BestBound and BestSol. Currently, your gap is 60% = (BestSol - BestBound) / BestSol

HiGHS will stop searching when the Gap is small (< 0.1%). You may also be happy with a larger gap, like 5%. You can change this by setting the mip_rel_gap parameter in HiGHS.

For your model, HiGHS has reached a point where it is no longer finding improving BestSol. This should give reasonable confidence that 0.5 is very nearly an optimal solution. However, HiGHS is struggling to improve the BestBound to prove that the BestSol is optimal.

If you don’t need a proof of optimality, you could interrupt the solve with CTRL+C and take the best solution. Or you could set a time limit (e.g., for 1 hr = 3600 sec, or even much shorter, like 60 sec).

If you do need a proof of optimality, you’ll need to improve your model, but this can be quite complicated if you don’t have a strong background in mixed-integer programming.

1 Like
7

how do you do this?

1 Like
8

JuMP.set_time_limit_sec(model, 60)

See JuMP · JuMP

1 Like
9

This is what I get

...
JuMP.set_time_limit_sec(model, 120)
...
# peso = Dict{String, Int}()
# for g in turni
#      if g in [t_festivi; "N"]
#         peso[g] = 2
#      else
#         peso[g] = 1
#      end
# end

peso = Dict{String, Float64}()
for g in turni
     if g in [t_festivi; "N"]
        peso[g] = 1
     else
        peso[g] = 0.5
     end
end

...


Symmetry detection completed in 0.0s
Found 31 full orbitope(s) acting on 486 columns

      2962     215        99   0.00%   0.1383315577    0.4999999996      72.33%     1095     17    508    183172    59.5s
 B    3049     252       129   0.00%   0.1383315577    0.4999999995      72.33%     1119     18    571    186656    60.2s
      3195     262       193   0.00%   0.1383315577    0.4999999995      72.33%     1155     19   1082    201234    65.3s
 L    3444     323       290   0.00%   0.1383315577    0.4999993007      72.33%     1027     10   1206    208151    71.6s
      3682     424       344   0.00%   0.1383315577    0.4999993007      72.33%      930     29   1367    250458    77.6s
      3899     492       423   0.00%   0.1383315577    0.4999993007      72.33%      938     16   1531    267557    82.6s
      4266     570       572   0.00%   0.1383315577    0.4999993007      72.33%      947     24   1897    288779    91.2s
      4461     606       654   0.00%   0.1383315577    0.4999993007      72.33%      729     16   2057    303380    97.6s
      4688     697       710   0.01%   0.1383315577    0.4999993007      72.33%      734     34   2319    330637   106.3s

Restarting search from the root node
Model after restart has 1018 rows, 955 cols (953 bin., 0 int., 0 impl., 2 cont., 0 dom.fix.), and 10482 nonzeros

      4863       0         0   0.00%   0.1383315577    0.4999993007      72.33%       28      0      0    343593   111.0s
      4863       0         0   0.00%   0.1384673982    0.4999993007      72.31%       28      7      8    348774   111.3s

Symmetry detection completed in 0.0s
Found 31 full orbitope(s) acting on 486 columns

      5847     221       378   0.00%   0.1384673982    0.4999993007      72.31%     1727     13    402    375026   116.4s
      6606     403       669   0.00%   0.1384673982    0.4999993007      72.31%     1600     27    958    399779   120.0s

Solving report
  Status            Time limit reached
  Primal bound      0.499999300654
  Dual bound        0.138467398201
  Gap               72.31% (tolerance: 0.01%)
  P-D integral      90.3413485309
  Solution status   feasible
                    0.499999300654 (objective)
                    0 (bound viol.)
                    2.98272517796e-12 (int. viol.)
                    0 (row viol.)
  Timing            120.01 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
                    13934 (separation)
                    135771 (heuristics)
ERROR: The model was not solved correctly. Here is the output of `solution_summary` to help debug why this happened:

solution_summary(; result = 1, verbose = false)
├ solver_name          : HiGHS
├ Termination
│ ├ termination_status : TIME_LIMIT
│ ├ result_count       : 1
│ ├ raw_status         : kHighsModelStatusTimeLimit
│ └ objective_bound    : 1.38467e-01
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : NO_SOLUTION
│ ├ objective_value      : 4.99999e-01
│ ├ dual_objective_value : NaN
│ └ relative_gap         : 7.23065e-01
└ Work counters
  ├ solve_time (sec)   : 1.20036e+02
  ├ simplex_iterations : 399779
  ├ barrier_iterations : -1
  └ node_count         : 6606

Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:35
 [2] #assert_is_solved_and_feasible#111
   @ C:\Users\sprmn\.julia\packages\JuMP\N7h14\src\optimizer_interface.jl:1008 [inlined]
 [3] assert_is_solved_and_feasible
   @ C:\Users\sprmn\.julia\packages\JuMP\N7h14\src\optimizer_interface.jl:1002 [inlined]
 [4] sol_jump(mese::Int64, y::Int64)
   @ Main c:\Users\sprmn\.julia\envs\v1_10\J_turni_SG_MWE.jl:236
 [5] top-level scope
   @ c:\Users\sprmn\.julia\envs\v1_10\J_turni_SG_MWE.jl:311

Instead if I use the version with integer weights, I quickly get a solution.

I haven’t thought about it very carefully, but at first glance I’d say that, since everything is linear, the two problems with proportional weights are equivalent.
So instead of weight1=1 and weight2=1.5,
I’ll use weight1=2 and weight2=3, which seems more digestible for HIGHS.

10

[3] assert_is_solved_and_feasible

Yes, you can’t use this function because HiGHS did not “solve” it to optimality. You should check if termination_status(model) == TIME_LIMIT and primal_status(model) == FEASIBLE_POINT.

but at first glance I’d say that, since everything is linear, the two problems with proportional weights are equivalent.

:+1: Yes, keeping data integer is usually helpful.

1 Like
11

With this line removed, here’s what I get.
The output tables are different (but “close”) to those obtained with the integer weights and the script that completed but not due to the time limit reached.


Solving report
  Status            Time limit reached
  Primal bound      0.499999300654
  Dual bound        0.138467398201
  Gap               72.31% (tolerance: 0.01%)
  P-D integral      90.2967188094
  Solution status   feasible
                    0.499999300654 (objective)
                    0 (bound viol.)
                    2.98272517796e-12 (int. viol.)
                    0 (row viol.)
  Timing            120.01 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
  Max sub-MIP depth 7
  Nodes             6717
  Repair LPs        0 (0 feasible; 0 iterations)
  LP iterations     401748 (total)
                    93517 (strong br.)
                    13941 (separation)
                    135771 (heuristics)
(3-dimensional DenseAxisArray{VariableRef,3,...} with index sets:
    Dimension 1, Base.OneTo(13)
    Dimension 2, Base.OneTo(31)
    Dimension 3, ["M1", "M2", "P", "N", "GF", "NF", "CUM1", "CUM2", "CLM", "CUMs", "CLMs"]
And data, a 13×31×11 Array{VariableRef, 3}:
[:, :, "M1"] =
...
1 Like
12

Sure. So the question now is: is that solution acceptable to you?

It’s probably good, but we can’t prove that there isn’t a solution that exists with a better objective value.

13

Very interesting. How do you change tolerance?

14

set_attribute(model, "mip_rel_gap", 0.05) (for 5%)

15

In this case the tolerance criterion intervened before the time limit






...


model = Model(HiGHS.Optimizer)
JuMP.set_time_limit_sec(model, 600)
set_attribute(model, "mip_rel_gap", 0.75)

...



Solving report
  Status            Optimal
  Primal bound      0.499999999996
  Dual bound        0.125
  Gap               75% (tolerance: 75%)
  P-D integral      10.7619481432
  Solution status   feasible
                    0.499999999996 (objective)
                    0 (bound viol.)
                    1.08695444656e-11 (int. viol.)
                    0 (row viol.)
  Timing            15.10 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
  Max sub-MIP depth 6
  Nodes             575
  Repair LPs        0 (0 feasible; 0 iterations)
  LP iterations     69878 (total)
                    485 (strong br.)
                    3060 (separation)
                    17260 (heuristics)
(3-dimensional DenseAxisArray{VariableRef,3,...} with index sets:
    Dimension 1, Base.OneTo(13)
    Dimension 2, Base.OneTo(31)
    Dimension 3, ["M1", "M2", "P", "N", "GF", "NF", "CUM1", "CUM2", "CLM", "CUMs", "CLMs"]
And data, a 13×31×11 Array{VariableRef, 3}:
[:, :, "M1"] =
 x[1,1,M1]   x[1,2,M1]   x[1,3,M1]   x[1,4,M1]   x[1,5,M1]   …  x[1,27,M1]   x[1,28,M1]   x[1,29,M1]   x[1,30,M1]   x[1,31,M1]
 x[2,1,M1]   x[2,2,M1]   x[2,3,M1]   x[2,4,M1]   x[2,5,M1]      x[2,27,M1]   x[2,28,M1]   x[2,29,M1]   x[2,30,M1]   x[2,31,M1]
 x[3,1,M1]   x[3,2,M1]   x[3,3,M1]   x[3,4,M1]   x[3,5,M1]
16

Yes, because you can see it terminates with OPTIMAL and Gap = 75%.

If I have one suggestion for you, is that you may at first focus on smaller scale systems (try a few toy examples) to build up your comprehension about the physical problem at hand, and the procedures of optimization.

I think currently you are involving too much staff
using XLSX, Dates , CSV, Tables, DataFrames, JuMP, HiGHS
that may distract you from focusing on more essential concepts. (e.g. the Gap, primal/dual bounds)

You can write a math program (e.g. using JuMP) from scratch at first. And after you fully understand your code, you may then add more decorations.

17

Yes, the problem is indeed complex. I can’t simplify the boundary conditions or the set of constraints without distorting the problem.
In this specific case, CSV and Tables don’t matter: they’re the remnants of previous attempts to address similar problems.
Other ancillary functions for calculating certain sets of days could also be removed or better consolidated.
But my goal, in this case, is purely practical. I need a solution, even an approximate one, that satisfies all the conditions mentioned.
In general, it’s nice to have a good understanding of how things work, and thanks to @odow, I’ve learned some things that, besides being theoretically interesting, I believe will be concretely useful in similar situations.

1 Like