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)