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

[Gagan_Meena]

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)

[odow]

See the tips and suggestions here:

• Solutions · JuMP
• Debugging · JuMP