# How to Find the constraints which are causing infeasibility in the model using compute_iis in Gurobi

How to compute the IIS form for the model shown below

``````using JuMP, Gurobi
include("Data.jl")

UnitCommitment_dt = UC_dt

function DUC_Clearing(UnitCommitment_dt)
# Indices
ngen = length(UC_dt.Gen_Constraints[:,1])
periods = length(UC_dt.Forecast[1,:])

# Sets
J = collect(1:ngen)
T = collect(1:periods)
zT = collect(1:periods+1)
# Parameters
# Generator Parameters
cแต = UC_dt.Gen_Constraints[:,1]
c  = UC_dt.Gen_Constraints[:,2]
Pฬฒ  = UC_dt.Gen_Constraints[:,3]
Pฬ  = UC_dt.Gen_Constraints[:,4]
Sแต = UC_dt.Gen_Constraints[:,5]
Sแดฐ = UC_dt.Gen_Constraints[:,6]
Rแต = UC_dt.Gen_Constraints[:,7]
Rแดฐ = UC_dt.Gen_Constraints[:,8]
Tแต = UC_dt.Gen_Constraints[:,9]
Tแดฐ = UC_dt.Gen_Constraints[:,10]

# Initial Init_Conditions
๐ฃโ = UC_dt.Init_Conditions[:,1]
๐โ = UC_dt.Init_Conditions[:,2]
U = UC_dt.Init_Conditions[:,3]
D = UC_dt.Init_Conditions[:,4]
L = min.(periods*ones(Real, 3), U)
F = min.(periods*ones(Real, 3), D)
๐งโโโ = [0; 0 ;0]

# Forecast
๐ท = UC_dt.Forecast[1,:]
๐ = UC_dt.Forecast[2,:]

# Model
DUC = Model(Gurobi.Optimizer)
set_attribute(DUC, "OutputFlag", 1)
set_attribute(DUC, "FeasibilityTol", 1e-8)
set_attribute(DUC, "MIPGap", 1e-8)
set_attribute(DUC, "MIPGapAbs", 0)

# Variables

@variables(DUC, begin
๐[J,T]          # pโฑผ(t): Active/real power produced by unit j in time period t (MW)
๐ฬ[J,T]          # ฬpโฑผ(t): Maximum available power in time period t from unit j (MW).
๐ฃ[J,T], Bin     # vโฑผ(t): On/off status in time period t of the unit j (1 if ON and 0 otherwise)
๐ฆ[J,T], Bin     # yโฑผ(t): Startup status in time period t of the unit j
๐ง[J,zT], Bin    # zโฑผ(t): Shutdown status in the time period t of the unit j
end)

@constraint(DUC, pow_balance[t=T], sum(๐[j,t] for j in J) == ๐ท[t])
@constraint(DUC, reserve[t=T], sum(๐ฬ[j,t] for j in J) โฅ ๐ท[t]+๐[t])

for t in T, j in J
if t!=1
@constraint(DUC, ๐ฃ[j,t-1]-๐ฃ[j,t]+๐ฆ[j,t]-๐ง[j,t]==0)
else
@constraint(DUC, ๐ฃโ[j]-๐ฃ[j,t]+๐ฆ[j,t]-๐ง[j,t]==0)
end
end
for t in T, j in J
if t!=1
@constraint(DUC, ๐[j,t]-๐[j,t-1] โค Rแต[j]*๐ฃ[j,t-1]+Sแต[j]*๐ฆ[j,t])
else
@constraint(DUC, ๐[j,t]-๐โ[j] โค Rแต[j]*๐ฃโ[j]+Sแต[j]*๐ฆ[j,t])
end
end
for t in T, j in J
if t!=1
@constraint(DUC, ๐[j,t-1]-๐[j,t] โค Rแดฐ[j]*๐ฃ[j,t]+Sแดฐ[j]*๐ง[j,t])
else
@constraint(DUC, ๐โ[j]-๐[j,t] โค Rแดฐ[j]*๐ฃ[j,t]+Sแดฐ[j]*๐ง[j,t])
end
end
@constraint(DUC, [j=J,t=T], Pฬฒ[j]*๐ฃ[j,t] โค ๐[j,t])
@constraint(DUC, [j=J,t=T], ๐[j,t] โค ๐ฬ[j,t])
@constraint(DUC, [j=J,t=T], ๐ฬ[j,t] โค Pฬ[j]*๐ฃ[j,t])
for t in T, j in J
if t!=1
@constraint(DUC, ๐ฬ[j,t] โค ๐[j,t-1] + Rแต[j]*๐ฃ[j,t-1]+Sแต[j]*๐ฆ[j,t])
else
@constraint(DUC, ๐ฬ[j,t] โค ๐โ[j] + Rแต[j]*๐ฃโ[j]+Sแต[j]*๐ฆ[j,t])
end
end
for t in T, j in J
#if t!=length(UC_dt.Forecast[1,:])
@constraint(DUC, ๐ฬ[j,t] โค Pฬ[j]*(๐ฃ[j,t]-๐ง[j,t+1]) + Sแดฐ[j]*๐ง[j,t+1])
#=
else
#@constraint(DUC, ๐ฬ[j,t] โค Pฬ[j]*(๐ฃ[j,t]-๐ง[j,1]) + Sแดฐ[j]*๐ง[j,1])
@constraint(DUC, ๐ฬ[j,t] โค Pฬ[j]*(๐ฃ[j,t]-๐ง[j,t+1]) + Sแดฐ[j]*๐ง[j,t+1])
end
=#

end
#constraint(DUC, ๐ฃ[3,1]==1)
for j in J
for t in L[j]+1:periods
k=t-Tแต[j]+1
ku = collect(k:t)
println(k)
println(ku)
if kโฅ1
@constraint(DUC, sum(๐ฆ[j,kuu] for kuu in ku) โค ๐ฃ[j,t])
end
end
end

for j in J
for t in F[j]+1:periods
k=t-Tแดฐ[j]+1
kd = collect(k:t)
if kโฅ1
@constraint(DUC, ๐ฃ[j,t] + sum(๐ง[j,kdd] for kdd in kd) โค 1)
end
end
end

for j in J
@assert(L[j]*F[j]==0)
for t in 1:L[j]
@constraint(DUC,๐ฃ[j,t]==1)
#@constraint(DUC,๐ฆ[j,t]==0)
#@constraint(DUC,๐ง[j,t]==0)
end
for t in 1:F[j]
@constraint(DUC,๐ฃ[j,t]==0)
#@constraint(DUC,๐ฆ[j,t]==0)
#@constraint(DUC,๐ง[j,t]==0)
end
end

@expressions(DUC, begin
gen_cost, sum(sum(c[j]*๐[j,t] for j in J) for t in T)
startup_cost, sum(sum(cแต[j]*๐ฆ[j,t] for j in J) for t in T)
end)
@objective(DUC, Min, gen_cost + startup_cost)
optimize!(DUC)

println(DUC)
println(value.(๐))
#println(dual.(pow_balance)) # can not find dual for MILP
println(dual_status(DUC))

# Finding dual for the linear version of the problem by fixing the discrete variables
undo = fix_discrete_variables(DUC)
println(" Here is the LP version of the model")
#println(DUC)

optimize!(DUC)
dual_status(DUC)
println(dual.(pow_balance))
println(dual.(reserve))
println(value.(๐ฬ))
println(value.(๐))
#println(value.(๐ง))

end

@time DUC_Clearing(UnitCommitment_dt)
``````

See the tips and suggestions here:

1 Like