CPLEX ignores a Initial feasible solution

Dear all,
I am solving a MILP with CPLEX 22.1. I have multiple objectives, and I minimize z1 first. I save the binary variables (called y[j,ℓ,t,τ]). Let Y0 the matriz with all y’s produced in the first stage.
I the second stage, I minimize z2, forcing z1 <= z1* + 0.1, by inserting Y0 as incumbent (only the binaries), to reduce the effort to solve the second problem. But, CPLEX ignores my solution (and it is feasible!).
When I fix these values at start value, I have a feasible problem.
I have used the following code to do this:

 for (j,ℓ,t,τ) in eachindex(y)
        set_start_value(y[j,ℓ,t,τ],Int(round(Y[j,ℓ,t,τ])))
 end

Anyone can help me please?

Hi @angeloaliano1,

Do you have the CPLEX log? How do you know that CPLEX ignored your solution?

Are there other discrete variables other than those binaries?

Could you try setting a primal start for all primal variables?

Hello @odow
Thanks for your response.
I will detail my problem.

  1. First phase: minimizing z1 (ok, works well)
λ = [0,0,1];
#----------------------------------------------------------------------------------
modelo = Model(CPLEX.Optimizer)
set_optimizer_attribute(modelo, "CPXPARAM_TimeLimit", 300.0) 
set_optimizer_attribute(modelo, "CPXPARAM_MIP_Tolerances_MIPGap", 0.001)
set_optimizer_attribute(modelo, "CPXPARAM_ScreenOutput", 1)

#Definindo o objetivo 1
@constraint(modelo, 
z_1 == 
sum(FC[ℓ,t]*y[j,ℓ,t,τ] for j in Je for ℓ in J_n_j[j] for t in 1:per-n[j]+1 for τ in t+n[j]-1:minimum([per,t+N[j]-1]) if γ[j] > 0) +
sum(RC[j,ℓ,t]*w[j,ℓ,t] for j in Je for ℓ in J_n_j[j] for t in T if γ[j] > 0) +
sum(OC[j,t]*(Q[j] - sum(w[j,ℓ,τ] for ℓ in J_n_j[j] for τ in 1:t)) for j in Je for t in T if γ[j] > 0) +
sum(OC[j,t]*Q[j] for j in Je for t in T if β[j] == 0) +
sum(OC[ℓ,t]*w[j,ℓ,τ] for ℓ in Jn  for t in 2:per for τ in 1:t-1 for j in J_e_l[ℓ] if γ[j] > 0) +
sum(UC[j,t]*u[j,t] for j in Je for t in T) +
sum(PC[i,r,t]*z1[i,r,j,m,t] for i in I for r in R for j in J for m in 1:2 for t in T if A[i,j] >= m) +
sum(PC[i,r,t]*z2[i,r,h,m,t] for i in I for r in R for h in H for m in 1:2 for t in T if A[i,h] >= m) + 
sum(MC[j,p,t]*x1[j,p,k,m,t] for j in J for p in P for k in K for m in 1:2 for t in T if A[j,k] >= m) + 
sum(MC[j,p,t]*x2[j,p,h,m,t] for j in J for p in P for h in H for m in 1:2 for t in T if A[j,h] >= m) +
sum(TC[i,j,m,t]*z1[i,r,j,m,t] for i in I for j in J for r in R for m in 1:2 for t in T if A[i,j] >= m) +
sum(TC[i,h,m,t]*z2[i,r,h,m,t] for i in I for h in H for r in R for m in 1:2 for t in T if A[i,h] >= m) +
sum(TC[h,hl,3,t]*z3[h,r,hl,t] for h in H for r in R for hl in H for t in T if A[h,hl] == 3) +
sum(TC[h,j,m,t]*z4[h,r,j,m,t] for h in H for r in R for j in J for m in 1:2 for t in T if A[h,j] >= m) +
sum(TC[j,k,m,t]*x1[j,p,k,m,t] for j in J for p in P for k in K for m in 1:2 for t in T if A[j,k] >= m) +
sum(TC[j,h,m,t]*x2[j,p,h,m,t] for j in J for p in P for h in H for m in 1:2 for t in T if A[j,h] >= m) +
sum(TC[h,hl,3,t]*x3[h,p,hl,t] for h in H for p in P for hl in H for t in T if A[h,hl] == 3) +
sum(TC[h,k,m,t]*x4[h,p,k,m,t] for h in H for p in P for k in K for m in 1:2 for t in T if A[h,k] >= m)                                                                              #Custo de operação do que sobrou
);
#Definindo o objetivo 2
@constraint(modelo,
z_2 == sum(e[i,j,m]*z1[i,r,j,m,t] for i in I for r in R for j in J for m in 1:2 for t in T if A[i,j] >= m) +                                #Emissões dos fornecedores aos serviços 
        sum(e[i,h,m]*z2[i,r,h,m,t] for i in I for r in R for h in H for m in 1:2 for t in T if A[i,h] >= m) +
        sum(e[h,hl,3]*z3[h,r,hl,t] for h in H for r in R for hl in H for t in T if A[h,hl] == 3) +
        sum(e[h,j,m]*z4[h,r,j,m,t] for h in H for r in R for j in J for m in 1:2 for t in T if A[h,j] >= m) +
        sum(e[j,k,m]*x1[j,p,k,m,t] for j in J for p in P for k in K for m in 1:2 for t in T if A[j,k] >= m) +                                #Emissões dos serviços aos clientes 
        sum(e[j,h,m]*x2[j,p,h,m,t] for j in J for p in P for h in H for m in 1:2 for t in T if A[j,h] >= m) +                                #Emissões de transporte das instalações até os portos
        sum(e[h,hl,3]*x3[h,p,hl,t] for h in H for p in P for hl in H for t in T if A[h,hl] == 3) +                                           #Emissões de transporte das entre os portos 
        sum(e[h,k,m]*x4[h,p,k,m,t] for h in H for p in P for k in K for m in 1:2 for t in T if A[h,k] >= m)                                  #Emissões de transporte dos portos até os clientes
);
#Definindo o objetivo 3
@constraint(modelo, 
    z_3 == -1*sum(w[j,ℓ,τ] for j in Je for ℓ in J_n_j[j] for t in 2:per for τ in 1:t if γ[j] > 0)
);
#A relocalização de uma instalação existente para exatamente um novo local 

@constraint(modelo,[j in Je; γ[j] > 0], 
    sum(y[j,ℓ,t,τ] for ℓ in J_n_j[j], t in 1:per-n[j]+1, τ in t+n[j]-1:minimum([per,t+N[j]-1])) == 1
);

#Estipulamos que cada nova facilidade reinstalada pode receber somente capacidade de uma única instalação existente
@constraint(modelo,con[ℓ in Jn], 
    sum(y[j,ℓ,t,τ] for j in J_e_l[ℓ], t in 1:per-n[j]+1, τ in t+n[j]-1:minimum([per,t+N[j]-1]) if γ[j] > 0) <= 1
);

@constraint(modelo, [j in Je, ℓ in J_n_j[j], t in T; γ[j] > 0], 
    w[j,ℓ,t] <= β[j]*Q[j]*sum(y[j,ℓ,τ,τl] for τ in maximum([1,t-N[j]+1]):minimum([per-n[j]+1,t]) for 
    τl in maximum([t,τ+n[j]-1]):minimum([per,τ+N[j]-1]))
);

#Garantem que a capacidade de transferência pode somente ocorrer durante o seu período de transferência
@constraint(modelo, [j in Je; γ[j] > 0], 
    sum(w[j,ℓ,t] for ℓ in J_n_j[j] for t in T) <= β[j]*Q[j]
);

...other constraints....


@constraint(modelo,[j in Je,t in T; γ[j] > 0],
    δ[j,t] <= sum(y[jl,ℓ,τ,τl] for jl in Je for ℓ in J_n_j[jl] for τ in maximum([1,t-N[jl]+1]):minimum([per-n[jl]+1,t]) for τl in maximum([t,τ+n[jl]-1]):minimum([per,τ+N[jl]-1]) if γ[jl] > 0 && jl != j)
)

@constraint(modelo,[j in Je,t in T; γ[j] > 0],
    δ[j,t] >= (1/(ηe-1))*sum(y[jl,ℓ,τ,τl] for jl in Je for ℓ in J_n_j[jl] for τ in maximum([1,t-N[jl]+1]):minimum([per-n[jl]+1,t]) for τl in maximum([t,τ+n[jl]-1]):minimum([per,τ+N[jl]-1]) if γ[jl] > 0 && jl != j)   
)

optimize!(modelo)
z3_value = value(z_3)
Y = value.(y)
  1. Second phase: fix z3 at optimal value, and the minimize z1 (works well too):
λ = [1,0,0]

modelo = Model(CPLEX.Optimizer)
set_optimizer_attribute(modelo, "CPXPARAM_TimeLimit", 300.0) 
set_optimizer_attribute(modelo, "CPXPARAM_MIP_Tolerances_MIPGap", 0.001)
set_optimizer_attribute(modelo, "CPXPARAM_ScreenOutput", 1)

#Objetivo 1: custos
@variable(modelo, z_1);
#Objetivo 2: emissões de CO2
@variable(modelo, z_2);
#Objetivo 3: impacto social
@variable(modelo, z_3);
#Definindo-se um objetivo ponderado 
@objective(modelo, Min, λ[1]*z_1 + λ[2]*z_2 + λ[3]*z_3); 
#@objective(modelo, Min, 1); 
#Definindo o objetivo 1

for (j,ℓ,t,τ) in eachindex(y)
    set_start_value(y[j,ℓ,t,τ],Int(round(Y[j,ℓ,t,τ])))
end

@constraint(modelo, z_3 <= z3_value + 10)

@constraint(modelo, 
z_1 == 
sum(FC[ℓ,t]*y[j,ℓ,t,τ] for j in Je for ℓ in J_n_j[j] for t in 1:per-n[j]+1 for τ in t+n[j]-1:minimum([per,t+N[j]-1]) if γ[j] > 0) +
sum(RC[j,ℓ,t]*w[j,ℓ,t] for j in Je for ℓ in J_n_j[j] for t in T if γ[j] > 0) +
sum(OC[j,t]*(Q[j] - sum(w[j,ℓ,τ] for ℓ in J_n_j[j] for τ in 1:t)) for j in Je for t in T if γ[j] > 0) +
sum(OC[j,t]*Q[j] for j in Je for t in T if β[j] == 0) +
sum(OC[ℓ,t]*w[j,ℓ,τ] for ℓ in Jn  for t in 2:per for τ in 1:t-1 for j in J_e_l[ℓ] if γ[j] > 0) +
sum(UC[j,t]*u[j,t] for j in Je for t in T) +
sum(PC[i,r,t]*z1[i,r,j,m,t] for i in I for r in R for j in J for m in 1:2 for t in T if A[i,j] >= m) +
sum(PC[i,r,t]*z2[i,r,h,m,t] for i in I for r in R for h in H for m in 1:2 for t in T if A[i,h] >= m) + 
sum(MC[j,p,t]*x1[j,p,k,m,t] for j in J for p in P for k in K for m in 1:2 for t in T if A[j,k] >= m) + 
sum(MC[j,p,t]*x2[j,p,h,m,t] for j in J for p in P for h in H for m in 1:2 for t in T if A[j,h] >= m) +
sum(TC[i,j,m,t]*z1[i,r,j,m,t] for i in I for j in J for r in R for m in 1:2 for t in T if A[i,j] >= m) +
sum(TC[i,h,m,t]*z2[i,r,h,m,t] for i in I for h in H for r in R for m in 1:2 for t in T if A[i,h] >= m) +
sum(TC[h,hl,3,t]*z3[h,r,hl,t] for h in H for r in R for hl in H for t in T if A[h,hl] == 3) +
sum(TC[h,j,m,t]*z4[h,r,j,m,t] for h in H for r in R for j in J for m in 1:2 for t in T if A[h,j] >= m) +
sum(TC[j,k,m,t]*x1[j,p,k,m,t] for j in J for p in P for k in K for m in 1:2 for t in T if A[j,k] >= m) +
sum(TC[j,h,m,t]*x2[j,p,h,m,t] for j in J for p in P for h in H for m in 1:2 for t in T if A[j,h] >= m) +
sum(TC[h,hl,3,t]*x3[h,p,hl,t] for h in H for p in P for hl in H for t in T if A[h,hl] == 3) +
sum(TC[h,k,m,t]*x4[h,p,k,m,t] for h in H for p in P for k in K for m in 1:2 for t in T if A[h,k] >= m)                                                                              #Custo de operação do que sobrou
);
#Definindo o objetivo 2
@constraint(modelo,
z_2 == sum(e[i,j,m]*z1[i,r,j,m,t] for i in I for r in R for j in J for m in 1:2 for t in T if A[i,j] >= m) +                                #Emissões dos fornecedores aos serviços 
        sum(e[i,h,m]*z2[i,r,h,m,t] for i in I for r in R for h in H for m in 1:2 for t in T if A[i,h] >= m) +
        sum(e[h,hl,3]*z3[h,r,hl,t] for h in H for r in R for hl in H for t in T if A[h,hl] == 3) +
        sum(e[h,j,m]*z4[h,r,j,m,t] for h in H for r in R for j in J for m in 1:2 for t in T if A[h,j] >= m) +
        sum(e[j,k,m]*x1[j,p,k,m,t] for j in J for p in P for k in K for m in 1:2 for t in T if A[j,k] >= m) +                                #Emissões dos serviços aos clientes 
        sum(e[j,h,m]*x2[j,p,h,m,t] for j in J for p in P for h in H for m in 1:2 for t in T if A[j,h] >= m) +                                #Emissões de transporte das instalações até os portos
        sum(e[h,hl,3]*x3[h,p,hl,t] for h in H for p in P for hl in H for t in T if A[h,hl] == 3) +                                           #Emissões de transporte das entre os portos 
        sum(e[h,k,m]*x4[h,p,k,m,t] for h in H for p in P for k in K for m in 1:2 for t in T if A[h,k] >= m)                                  #Emissões de transporte dos portos até os clientes
);
#Definindo o objetivo 3
@constraint(modelo, 
    z_3 == -1*sum(w[j,ℓ,τ] for j in Je for ℓ in J_n_j[j] for t in 2:per for τ in 1:t if γ[j] > 0)
);
#A relocalização de uma instalação existente para exatamente um novo local 

@constraint(modelo,[j in Je; γ[j] > 0], 
    sum(y[j,ℓ,t,τ] for ℓ in J_n_j[j], t in 1:per-n[j]+1, τ in t+n[j]-1:minimum([per,t+N[j]-1])) == 1
);

#Estipulamos que cada nova facilidade reinstalada pode receber somente capacidade de uma única instalação existente
@constraint(modelo,con[ℓ in Jn], 
    sum(y[j,ℓ,t,τ] for j in J_e_l[ℓ], t in 1:per-n[j]+1, τ in t+n[j]-1:minimum([per,t+N[j]-1]) if γ[j] > 0) <= 1
);

@constraint(modelo, [j in Je, ℓ in J_n_j[j], t in T; γ[j] > 0], 
    w[j,ℓ,t] <= β[j]*Q[j]*sum(y[j,ℓ,τ,τl] for τ in maximum([1,t-N[j]+1]):minimum([per-n[j]+1,t]) for 
    τl in maximum([t,τ+n[j]-1]):minimum([per,τ+N[j]-1]))
);

#Garantem que a capacidade de transferência pode somente ocorrer durante o seu período de transferência
@constraint(modelo, [j in Je; γ[j] > 0], 
    sum(w[j,ℓ,t] for ℓ in J_n_j[j] for t in T) <= β[j]*Q[j]
);

@constraint(modelo, [j in Je; γ[j] > 0], 
    sum(w[j,ℓ,t] for ℓ in J_n_j[j] for t in T) >= γ[j]*Q[j]
);

@constraint(modelo,[j in Je,ℓ in J_n_j[j], t in T; γ[j] > 0],
    w[j,ℓ,t] >= π[j]*Q[j]*sum(y[j,ℓ,τ,τl] for τ in maximum([1,t-N[j]+1]):minimum([per-n[j]+1,t]) for 
    τl in maximum([t,τ+n[j]-1]):minimum([per,τ+N[j]-1]))
)

#Capacidade do fornecedor
@constraint(modelo, [i in I, r in R, t in T], 
    sum(z1[i,r,j,m,t] for j in J for m in 1:2 if A[i,j] >= m) + 
    sum(z2[i,r,h,m,t] for h in H for m in 1:2 if A[i,h] >= m) <= Q_barra[i,r,t]
);

#Previnem que facilidades existentes excedem suas capacidades se estão envolvidas ou não na relocalização
@constraint(modelo, [j in Je, t in T; β[j] == 0], 
    sum(θ[p]*x1[j,p,k,m,t] for p in P for k in K for m in 1:2 if A[j,k] >= m) + 
    sum(θ[p]*x2[j,p,h,m,t] for p in P for h in H for m in 1:2 if A[j,h] >= m) <= Q[j] + u[j,t]
);

@constraint(modelo, [j in Je, t in T; γ[j] > 0], 
    sum(θ[p]*x1[j,p,k,m,t] for p in P for k in K for m in 1:2 if A[j,k] >= m) + 
    sum(θ[p]*x2[j,p,h,m,t] for p in P for h in H for m in 1:2 if A[j,h] >= m) <= 
    Q[j] + u[j,t] - sum(w[j,ℓ,τ] for ℓ in J_n_j[j] for τ in 1:t)
);

@constraint(modelo,[j in Jn,p in P,k in K,m in [1,2],t in [1]; A[j,k] >= m],
    x1[j,p,k,m,t] == 0 
)

@constraint(modelo,[j in Jn,p in P,h in H,m in [1,2],t in [1]; A[j,h] >= m],
    x2[j,p,h,m,t] == 0 
)

@constraint(modelo,
    sum(w[j,ℓ,per] for j in Je for ℓ in J_n_j[j] if β[j] > 0) == 0
)

#Condições similares às duas anteriores
@constraint(modelo, [ℓ in Jn, t in 2:per], 
    sum(θ[p]*x1[ℓ,p,k,m,t] for p in P for k in K for m in 1:2 if A[ℓ,k] >= m) + 
    sum(θ[p]*x2[ℓ,p,h,m,t] for p in P for h in H for m in 1:2 if A[ℓ,h] >= m) <= 
    sum(w[j,ℓ,τ] for j in J_e_l[ℓ] for τ in 1:t-1 if γ[j] > 0)
);

#Conservação de fluxo entre os fornecedores e as instalações
@constraint(modelo, [j in J, r in R, t in T],
    sum(z1[i,r,j,m,t] for i in I for m in 1:2 if A[i,j] >= m) + 
    sum(z4[h,r,j,m,t] for h in H for m in 1:2 if A[h,j] >= m) == 
    sum(a[r,p]*x1[j,p,k,m,t] for p in P for k in K for m in 1:2 if A[j,k] >= m) + 
    sum(a[r,p]*x2[j,p,h,m,t] for p in P for h in H for m in 1:2 if A[j,h] >= m)
);

#Conservação do fluxo das matérias primas nos portos (1)
@constraint(modelo, [r in R, h in H, t in T], 
    sum(z2[i,r,h,m,t] for i in I for m in 1:2 if A[i,h] >= m) == 
    sum(z3[h,r,hl,t] for hl in H if A[h,hl] == 3) 
);

#Conservação do fluxo das matérias primas nos portos (2)
@constraint(modelo, [r in R, h in H, t in T], 
    sum(z3[hl,r,h,t] for hl in H if A[hl,h] == 3) == 
    sum(z4[h,r,j,m,t] for j in J for m in 1:2 if A[h,j] >= m)
);

#Conservação do fluxo dos produtos nos portos (1) 
@constraint(modelo, [p in P, h in H, t in T], 
    sum(x2[j,p,h,m,t] for j in J for m in 1:2 if A[j,h] >= m) == 
    sum(x3[h,p,hl,t] for hl in H if A[h,hl] == 3) 
);
#Conservação do fluxo dos produtos nos portos (2)
@constraint(modelo, [p in P, h in H, t in T], 
    sum(x3[hl,p,h,t] for hl in H if A[hl,h] == 3) == 
    sum(x4[h,p,k,m,t] for k in K for m in 1:2 if A[h,k] >= m)
);

#Atendimento da demanda dos clientes
@constraint(modelo, [p in P, k in K, t in T], 
    sum(x1[j,p,k,m,t] for j in J for m in 1:2 if A[j,k] >= m) + 
    sum(x4[h,p,k,m,t] for h in H for m in 1:2 if A[h,k] >=m) == D[k,p,t] 
);

@constraint(modelo,[t in T],
    sum(u[j,t] for j in Je) <= sum(w[j,ℓ,t] for j in Je for ℓ in J_n_j[j] if γ[j] > 0)
)

@constraint(modelo,[j in Je, t in T; γ[j] > 0],
    u[j,t] <= Qa[j]*(1 - sum(y[j,ℓ,τ,τl] for ℓ in J_n_j[j] for τ in 1:minimum([per-n[j]+1,t]) for τl in τ+n[j]-1:minimum([per,τ+N[j]-1])))    
)

@constraint(modelo,[j in Je, t in T],
    u[j,t] <= Qa[j]*δ[j,t]
)

@constraint(modelo,[j in Je,t in T; β[j] == 0],
    δ[j,t] <= sum(y[jl,ℓ,τ,τl] for jl in Je for ℓ in J_n_j[jl] for τ in maximum([1,t-N[jl]+1]):minimum([per-n[jl]+1,t]) for τl in maximum([t,τ+n[jl]-1]):minimum([per,τ+N[jl]-1]) if γ[jl] > 0)
)

@constraint(modelo,[j in Je,t in T; β[j] == 0],
    δ[j,t] >= (1/ηe)*sum(y[jl,ℓ,τ,τl] for jl in Je for ℓ in J_n_j[jl] for τ in maximum([1,t-N[jl]+1]):minimum([per-n[jl]+1,t]) for τl in maximum([t,τ+n[jl]-1]):minimum([per,τ+N[jl]-1]) if γ[jl] > 0)
)

@constraint(modelo,[j in Je,t in T; γ[j] > 0],
    δ[j,t] <= sum(y[jl,ℓ,τ,τl] for jl in Je for ℓ in J_n_j[jl] for τ in maximum([1,t-N[jl]+1]):minimum([per-n[jl]+1,t]) for τl in maximum([t,τ+n[jl]-1]):minimum([per,τ+N[jl]-1]) if γ[jl] > 0 && jl != j)
)

@constraint(modelo,[j in Je,t in T; γ[j] > 0],
    δ[j,t] >= (1/(ηe-1))*sum(y[jl,ℓ,τ,τl] for jl in Je for ℓ in J_n_j[jl] for τ in maximum([1,t-N[jl]+1]):minimum([per-n[jl]+1,t]) for τl in maximum([t,τ+n[jl]-1]):minimum([per,τ+N[jl]-1]) if γ[jl] > 0 && jl != j)   
)

optimize!(modelo)
z3_value = value(z_3)
z1_value = value(z_1)
Y = value.(y)

cplex report:

julia> 
Version identifier: 22.1.1.0 | 2022-11-27 | 9160aff4d
CPXPARAM_TimeLimit                               300
CPXPARAM_MIP_Tolerances_MIPGap                   0.001
Tried aggregator 3 times.
MIP Presolve eliminated 1015 rows and 11163 columns.
MIP Presolve modified 304 coefficients.
Aggregator did 36 substitutions.
Reduced MIP has 4279 rows, 26412 columns, and 161193 nonzeros.
Reduced MIP has 638 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.14 sec. (178.82 ticks)
Found incumbent of value -7.7372392e+08 after 0.23 sec. (297.78 ticks)
Tried aggregator 1 time.
Detecting symmetries...
MIP Presolve eliminated 1508 rows and 15064 columns.
Reduced MIP has 2771 rows, 11348 columns, and 122093 nonzeros.
Reduced MIP has 638 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.17 sec. (225.88 ticks)
Probing time = 0.09 sec. (41.12 ticks)
Cover probing fixed 8 vars, tightened 2733 bounds.
Clique table members: 49927.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 12 threads.
Root relaxation solution time = 0.06 sec. (110.90 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

*     0+    0                      -7.73724e+08  -5.87938e+10              --- 
      0     0  -1.19016e+09     9  -7.73724e+08  -1.19016e+09     1899   53.82%
*     0+    0                      -1.19016e+09  -1.19016e+09             0.00%
      0     0        cutoff        -1.19016e+09  -1.19016e+09     1899    0.00%
Elapsed time = 0.67 sec. (777.02 ticks, tree = 0.01 MB, solutions = 2)

Root node processing (before b&c):
  Real time             =    0.67 sec. (778.73 ticks)
Parallel b&c, 12 threads:
  Real time             =    0.00 sec. (0.00 ticks)
  Sync time (average)   =    0.00 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =    0.67 sec. (778.73 ticks)
  1. Phase 3: fix z3 and z1 at the optimal values (from previous problems) providing the binary variables y’s as incumbent (works bad, because we don’t have any integer solution). We minimize z2:
λ = [0,1,0]

modelo = Model(CPLEX.Optimizer)
set_optimizer_attribute(modelo, "CPXPARAM_TimeLimit", 300.0) 
set_optimizer_attribute(modelo, "CPXPARAM_MIP_Tolerances_MIPGap", 0.001)
set_optimizer_attribute(modelo, "CPXPARAM_ScreenOutput", 1)

#Variável binária y[j,ℓ,t,τ] = 1 se uma facilidade existente j é reinstalada no local l começando no períoodo t e terminando no período τ e 0, caso contrário
@variable(modelo, y[j in Je, ℓ in J_n_j[j], t in 1:per-n[j]+1, τ in t+n[j]-1:minimum([per,t+N[j]-1]); γ[j] > 0], Bin);
#Quantidade de capacidade transferida do local j para o novo local ℓ no período t
@variable(modelo, w[j in Je,ℓ in J_n_j[j],t in T] >=0);
#Quantidade de matéria prima r transportada do fornecedor i para a facilidade j no período t usando o modal m=1,2 (rodo ou ferro)
@variable(modelo, z1[i in I,r in R,j in J,m in [1,2],t in T; A[i,j] >= m] >=0);
#Quantidade de matéria prima r transportada do fornecedor i para o porto h no período t usando o modal m=1,2 (rodo ou ferro)
@variable(modelo, z2[i in I,r in R,h in H,m in [1,2],t in T; A[i,h] >= m] >=0);
#Quantidade de matéria prima r transportada do porto h para o porto hl no período t (hidro apenas)
@variable(modelo, z3[h in H,r in R,hl in H,t in T; A[h,hl] == 3] >=0);
#Quantidade de matéria prima r transportada do porto h para a facilidade j no período t usando o modal m=1,2 (rodo ou ferro)
@variable(modelo, z4[h in H,r in R,j in J,m in [1,2],t in T; A[h,j] >= m] >=0);
#Quantidade de produto p transportada do serviço j para o cliente k no período t pelo modal m=1,2 (rodo ou ferro)
@variable(modelo, x1[j in J,p in P,k in K,m in [1,2],t in T; A[j,k] >= m] >=0);
#Quantidade de produto p transportada do serviço j para o porto h no período t pelo modal m=1,2 (rodo ou ferro)
@variable(modelo, x2[j in J,p in P,h in H,m in [1,2],t in T; A[j,h] >= m] >=0);
#Quantidade de produto p transportada do porto h para o porto hl no período t (hidro apenas)
@variable(modelo, x3[h in H,p in P,hl in H,t in T; A[h,hl] == 3] >=0);
#Quantidade de produto p transportada do porto h para o cliente k com o modal m no período t (rodo e ferro)
@variable(modelo, x4[h in H,p in P,k in K,m in [1,2],t in T; A[h,k] >= m] >=0);
#Binária que indica se a facilidade j é elegível para usar capacidade adicional no período t
@variable(modelo, δ[j in Je,t in T], Bin)
#Quantidade dee capacidade adicional usada pela facilidade j no período t 
@variable(modelo, u[j in Je,t in T] >=0)
#Objetivo 1: custos
@variable(modelo, z_1);
#Objetivo 2: emissões de CO2
@variable(modelo, z_2);
#Objetivo 3: impacto social
@variable(modelo, z_3);
#Definindo-se um objetivo ponderado 
@objective(modelo, Min, λ[1]*z_1 + λ[2]*z_2 + λ[3]*z_3); 
#@objective(modelo, Min, 1); 
#Definindo o objetivo 1

for (j,ℓ,t,τ) in eachindex(y)
    set_start_value(y[j,ℓ,t,τ],Int(round(Y[j,ℓ,t,τ])))
end

@constraint(modelo, z_1 <= z1_value + 10)
@constraint(modelo, z_3 <= z3_value + 10)

@constraint(modelo, 
z_1 == 
sum(FC[ℓ,t]*y[j,ℓ,t,τ] for j in Je for ℓ in J_n_j[j] for t in 1:per-n[j]+1 for τ in t+n[j]-1:minimum([per,t+N[j]-1]) if γ[j] > 0) +
sum(RC[j,ℓ,t]*w[j,ℓ,t] for j in Je for ℓ in J_n_j[j] for t in T if γ[j] > 0) +
sum(OC[j,t]*(Q[j] - sum(w[j,ℓ,τ] for ℓ in J_n_j[j] for τ in 1:t)) for j in Je for t in T if γ[j] > 0) +
sum(OC[j,t]*Q[j] for j in Je for t in T if β[j] == 0) +
sum(OC[ℓ,t]*w[j,ℓ,τ] for ℓ in Jn  for t in 2:per for τ in 1:t-1 for j in J_e_l[ℓ] if γ[j] > 0) +
sum(UC[j,t]*u[j,t] for j in Je for t in T) +
sum(PC[i,r,t]*z1[i,r,j,m,t] for i in I for r in R for j in J for m in 1:2 for t in T if A[i,j] >= m) +
sum(PC[i,r,t]*z2[i,r,h,m,t] for i in I for r in R for h in H for m in 1:2 for t in T if A[i,h] >= m) + 
sum(MC[j,p,t]*x1[j,p,k,m,t] for j in J for p in P for k in K for m in 1:2 for t in T if A[j,k] >= m) + 
sum(MC[j,p,t]*x2[j,p,h,m,t] for j in J for p in P for h in H for m in 1:2 for t in T if A[j,h] >= m) +
sum(TC[i,j,m,t]*z1[i,r,j,m,t] for i in I for j in J for r in R for m in 1:2 for t in T if A[i,j] >= m) +
sum(TC[i,h,m,t]*z2[i,r,h,m,t] for i in I for h in H for r in R for m in 1:2 for t in T if A[i,h] >= m) +
sum(TC[h,hl,3,t]*z3[h,r,hl,t] for h in H for r in R for hl in H for t in T if A[h,hl] == 3) +
sum(TC[h,j,m,t]*z4[h,r,j,m,t] for h in H for r in R for j in J for m in 1:2 for t in T if A[h,j] >= m) +
sum(TC[j,k,m,t]*x1[j,p,k,m,t] for j in J for p in P for k in K for m in 1:2 for t in T if A[j,k] >= m) +
sum(TC[j,h,m,t]*x2[j,p,h,m,t] for j in J for p in P for h in H for m in 1:2 for t in T if A[j,h] >= m) +
sum(TC[h,hl,3,t]*x3[h,p,hl,t] for h in H for p in P for hl in H for t in T if A[h,hl] == 3) +
sum(TC[h,k,m,t]*x4[h,p,k,m,t] for h in H for p in P for k in K for m in 1:2 for t in T if A[h,k] >= m)                                                                              #Custo de operação do que sobrou
);
#Definindo o objetivo 2
@constraint(modelo,
z_2 == sum(e[i,j,m]*z1[i,r,j,m,t] for i in I for r in R for j in J for m in 1:2 for t in T if A[i,j] >= m) +                                #Emissões dos fornecedores aos serviços 
        sum(e[i,h,m]*z2[i,r,h,m,t] for i in I for r in R for h in H for m in 1:2 for t in T if A[i,h] >= m) +
        sum(e[h,hl,3]*z3[h,r,hl,t] for h in H for r in R for hl in H for t in T if A[h,hl] == 3) +
        sum(e[h,j,m]*z4[h,r,j,m,t] for h in H for r in R for j in J for m in 1:2 for t in T if A[h,j] >= m) +
        sum(e[j,k,m]*x1[j,p,k,m,t] for j in J for p in P for k in K for m in 1:2 for t in T if A[j,k] >= m) +                                #Emissões dos serviços aos clientes 
        sum(e[j,h,m]*x2[j,p,h,m,t] for j in J for p in P for h in H for m in 1:2 for t in T if A[j,h] >= m) +                                #Emissões de transporte das instalações até os portos
        sum(e[h,hl,3]*x3[h,p,hl,t] for h in H for p in P for hl in H for t in T if A[h,hl] == 3) +                                           #Emissões de transporte das entre os portos 
        sum(e[h,k,m]*x4[h,p,k,m,t] for h in H for p in P for k in K for m in 1:2 for t in T if A[h,k] >= m)                                  #Emissões de transporte dos portos até os clientes
);
#Definindo o objetivo 3
@constraint(modelo, 
    z_3 == -1*sum(w[j,ℓ,τ] for j in Je for ℓ in J_n_j[j] for t in 2:per for τ in 1:t if γ[j] > 0)
);
#A relocalização de uma instalação existente para exatamente um novo local 

@constraint(modelo,[j in Je; γ[j] > 0], 
    sum(y[j,ℓ,t,τ] for ℓ in J_n_j[j], t in 1:per-n[j]+1, τ in t+n[j]-1:minimum([per,t+N[j]-1])) == 1
);

#Estipulamos que cada nova facilidade reinstalada pode receber somente capacidade de uma única instalação existente
@constraint(modelo,con[ℓ in Jn], 
    sum(y[j,ℓ,t,τ] for j in J_e_l[ℓ], t in 1:per-n[j]+1, τ in t+n[j]-1:minimum([per,t+N[j]-1]) if γ[j] > 0) <= 1
);

@constraint(modelo, [j in Je, ℓ in J_n_j[j], t in T; γ[j] > 0], 
    w[j,ℓ,t] <= β[j]*Q[j]*sum(y[j,ℓ,τ,τl] for τ in maximum([1,t-N[j]+1]):minimum([per-n[j]+1,t]) for 
    τl in maximum([t,τ+n[j]-1]):minimum([per,τ+N[j]-1]))
);

#Garantem que a capacidade de transferência pode somente ocorrer durante o seu período de transferência
@constraint(modelo, [j in Je; γ[j] > 0], 
    sum(w[j,ℓ,t] for ℓ in J_n_j[j] for t in T) <= β[j]*Q[j]
);

@constraint(modelo, [j in Je; γ[j] > 0], 
    sum(w[j,ℓ,t] for ℓ in J_n_j[j] for t in T) >= γ[j]*Q[j]
);

@constraint(modelo,[j in Je,ℓ in J_n_j[j], t in T; γ[j] > 0],
    w[j,ℓ,t] >= π[j]*Q[j]*sum(y[j,ℓ,τ,τl] for τ in maximum([1,t-N[j]+1]):minimum([per-n[j]+1,t]) for 
    τl in maximum([t,τ+n[j]-1]):minimum([per,τ+N[j]-1]))
)

#Capacidade do fornecedor
@constraint(modelo, [i in I, r in R, t in T], 
    sum(z1[i,r,j,m,t] for j in J for m in 1:2 if A[i,j] >= m) + 
    sum(z2[i,r,h,m,t] for h in H for m in 1:2 if A[i,h] >= m) <= Q_barra[i,r,t]
);

#Previnem que facilidades existentes excedem suas capacidades se estão envolvidas ou não na relocalização
@constraint(modelo, [j in Je, t in T; β[j] == 0], 
    sum(θ[p]*x1[j,p,k,m,t] for p in P for k in K for m in 1:2 if A[j,k] >= m) + 
    sum(θ[p]*x2[j,p,h,m,t] for p in P for h in H for m in 1:2 if A[j,h] >= m) <= Q[j] + u[j,t]
);

@constraint(modelo, [j in Je, t in T; γ[j] > 0], 
    sum(θ[p]*x1[j,p,k,m,t] for p in P for k in K for m in 1:2 if A[j,k] >= m) + 
    sum(θ[p]*x2[j,p,h,m,t] for p in P for h in H for m in 1:2 if A[j,h] >= m) <= 
    Q[j] + u[j,t] - sum(w[j,ℓ,τ] for ℓ in J_n_j[j] for τ in 1:t)
);

@constraint(modelo,[j in Jn,p in P,k in K,m in [1,2],t in [1]; A[j,k] >= m],
    x1[j,p,k,m,t] == 0 
)

@constraint(modelo,[j in Jn,p in P,h in H,m in [1,2],t in [1]; A[j,h] >= m],
    x2[j,p,h,m,t] == 0 
)

@constraint(modelo,
    sum(w[j,ℓ,per] for j in Je for ℓ in J_n_j[j] if β[j] > 0) == 0
)

#Condições similares às duas anteriores
@constraint(modelo, [ℓ in Jn, t in 2:per], 
    sum(θ[p]*x1[ℓ,p,k,m,t] for p in P for k in K for m in 1:2 if A[ℓ,k] >= m) + 
    sum(θ[p]*x2[ℓ,p,h,m,t] for p in P for h in H for m in 1:2 if A[ℓ,h] >= m) <= 
    sum(w[j,ℓ,τ] for j in J_e_l[ℓ] for τ in 1:t-1 if γ[j] > 0)
);

...other constraints....

@constraint(modelo,[t in T],
    sum(u[j,t] for j in Je) <= sum(w[j,ℓ,t] for j in Je for ℓ in J_n_j[j] if γ[j] > 0)
)

@constraint(modelo,[j in Je, t in T; γ[j] > 0],
    u[j,t] <= Qa[j]*(1 - sum(y[j,ℓ,τ,τl] for ℓ in J_n_j[j] for τ in 1:minimum([per-n[j]+1,t]) for τl in τ+n[j]-1:minimum([per,τ+N[j]-1])))    
)

@constraint(modelo,[j in Je, t in T],
    u[j,t] <= Qa[j]*δ[j,t]
)

@constraint(modelo,[j in Je,t in T; β[j] == 0],
    δ[j,t] <= sum(y[jl,ℓ,τ,τl] for jl in Je for ℓ in J_n_j[jl] for τ in maximum([1,t-N[jl]+1]):minimum([per-n[jl]+1,t]) for τl in maximum([t,τ+n[jl]-1]):minimum([per,τ+N[jl]-1]) if γ[jl] > 0)
)

@constraint(modelo,[j in Je,t in T; β[j] == 0],
    δ[j,t] >= (1/ηe)*sum(y[jl,ℓ,τ,τl] for jl in Je for ℓ in J_n_j[jl] for τ in maximum([1,t-N[jl]+1]):minimum([per-n[jl]+1,t]) for τl in maximum([t,τ+n[jl]-1]):minimum([per,τ+N[jl]-1]) if γ[jl] > 0)
)

@constraint(modelo,[j in Je,t in T; γ[j] > 0],
    δ[j,t] <= sum(y[jl,ℓ,τ,τl] for jl in Je for ℓ in J_n_j[jl] for τ in maximum([1,t-N[jl]+1]):minimum([per-n[jl]+1,t]) for τl in maximum([t,τ+n[jl]-1]):minimum([per,τ+N[jl]-1]) if γ[jl] > 0 && jl != j)
)

@constraint(modelo,[j in Je,t in T; γ[j] > 0],
    δ[j,t] >= (1/(ηe-1))*sum(y[jl,ℓ,τ,τl] for jl in Je for ℓ in J_n_j[jl] for τ in maximum([1,t-N[jl]+1]):minimum([per-n[jl]+1,t]) for τl in maximum([t,τ+n[jl]-1]):minimum([per,τ+N[jl]-1]) if γ[jl] > 0 && jl != j)   
)

optimize!(modelo)

cplex report:

Version identifier: 22.1.1.0 | 2022-11-27 | 9160aff4d
CPXPARAM_TimeLimit                               300
CPXPARAM_MIP_Tolerances_MIPGap                   0.001
Warning:  No solution found from 1 MIP starts.
Retaining values of one MIP start for possible repair.
Tried aggregator 2 times.
MIP Presolve eliminated 1015 rows and 3105 columns.
MIP Presolve modified 55964 coefficients.
Aggregator did 24 substitutions.
Reduced MIP has 4293 rows, 34482 columns, and 215410 nonzeros.
Reduced MIP has 638 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 3.41 sec. (2627.09 ticks)
Tried aggregator 1 time.
Detecting symmetries...
Reduced MIP has 4293 rows, 34482 columns, and 215410 nonzeros.
Reduced MIP has 638 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.44 sec. (137.75 ticks)
Probing time = 1.08 sec. (63.87 ticks)
Cover probing fixed 8 vars, tightened 846 bounds.
Clique table members: 49927.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 12 threads.
Root relaxation solution time = 11.66 sec. (2277.85 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

      0     0   1.00839e+10    11                 1.00839e+10    38687
      0     0   1.02988e+10    10                Flowcuts: 14    39443
      0     0   1.03667e+10    12                    Cuts: 17    40350
      0     0   1.05099e+10    11                Flowcuts: 17    41321
      0     0   1.06284e+10    10                Flowcuts: 10    41988
      0     0   1.06575e+10    13                Flowcuts: 13    42442
Detecting symmetries...
      0     0   1.06584e+10    11                Flowcuts: 24    42461
      0     0   1.06589e+10    11                 Flowcuts: 3    42486
      0     0   1.06589e+10    11                 Flowcuts: 4    42496
Repair heuristic found nothing.
Detecting symmetries...
      0     2   1.06589e+10    11                 1.06589e+10    42496
Elapsed time = 72.26 sec. (31388.98 ticks, tree = 0.02 MB, solutions = 0)
      2     4   1.06767e+10    10                 1.06590e+10    42983
      3     1    infeasible                       1.06590e+10    43707
      4     2    infeasible                       1.06590e+10    44104
      5     3    infeasible                       1.06590e+10    44641
      6     5   1.08788e+10     7                 1.06590e+10    57191
      7     1    infeasible                       1.09457e+10    62221
      8     1    infeasible                       1.10190e+10    62502

Flow cuts applied:  52

Root node processing (before b&c):
  Real time             =   71.98 sec. (30983.80 ticks)
Parallel b&c, 12 threads:
  Real time             =   24.27 sec. (17962.71 ticks)
  Sync time (average)   =   21.13 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =   96.25 sec. (48946.51 ticks)

The models are the same (I don’t have space to insert the three complete codes - but the models are almost the same), I only add the constraints imposing an upper bound the the objective previously optimized… Why problem 3 are infeasible? Could you help me please?

I can’t run your codes because I don’t have the data etc, so I can’t explain why the third problem is infeasible. The answer is probably that you have some typo that differs between the models.

Are you trying to solve a lexicographic multi objective problem?

Instead of writing out separate models, you can write a loop:

using JuMP, CPLEX
model = Model(CPLEX.Optimizer)
@variable(model, z[1:3])
# other model ...
for i in 1:2
    @objective(model, Min, z[i])
    optimize!(model)
    assert_is_solved_and_feasible(model)
    z_value = value(z[i])
    @constraint(model, z[i] <= z_value + 10)
end
@objective(model, Min, z[3])
optimize!(model)
assert_is_solved_and_feasible(model)

Have you seen MultiObjectiveAlgorithms.jl? You can alternatively do something like:

using JuMP, CPLEX
import MultiObjectiveAlgorithms as MOA
model = Model(() -> MOA.Optimizer(CPLEX.Optimizer))
@variable(model, z[1:3])
@objective(model, Min, z)
set_attribute(model, MOA.Algorithm(), MOA.Lexicographic())
set_attribute(model, MOA.LexicographicAllPermutations(), false)
set_attribute(model, MOA.ObjectiveRelativeTolerance(), 0.01) # obj allowed to increase by 1%
optimize!(model)
assert_is_solved_and_feasible(model)
z_opt = objective_value(model)

See also Simple multi-objective examples · JuMP