Hello @odow
Thanks for your response.
I will detail my problem.
- 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)
- 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)
- 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?