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)
# 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,T], 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])
end
end
for j in J
for t in L[j]+1:periods
k=t-Tᵁ[j]+1
kd = collect(k:t)
if k!=0
@constraint(DUC, sum(𝑦[j,kdd] for kdd in kd) ≤ 𝑣[j,t])
end
print(k)
print(t)
println("")
end
end
for j in J
for t in F[j]+1:periods
k=t-Tᴰ[j]+1
kd = collect(k:t)
if k!=0
@constraint(DUC, 𝑣[j,t] + sum(𝑧[j,kdd] for kdd in kd) ≤ 1)
end
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))
end
@time DUC_Clearing(UnitCommitment_dt)
The above code is showing the following error
Gurobi Error 10005: Unable to retrieve attribute 'Pi'
Stacktrace:
[1] _check_ret
@ C:\Users\Admin\.julia\packages\Gurobi\EKa6j\src\MOI_wrapper\MOI_wrapper.jl:375 [inlined]
[2] get(model::Gurobi.Optimizer, attr::MathOptInterface.ConstraintDual, c::MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}})
@ Gurobi C:\Users\Admin\.julia\packages\Gurobi\EKa6j\src\MOI_wrapper\MOI_wrapper.jl:3086
[3] get
@ C:\Users\Admin\.julia\packages\MathOptInterface\wrorD\src\Bridges\bridge_optimizer.jl:1424 [inlined]
[4] get(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Gurobi.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, attr::MathOptInterface.ConstraintDual, index::MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}})
@ MathOptInterface.Utilities C:\Users\Admin\.julia\packages\MathOptInterface\wrorD\src\Utilities\cachingoptimizer.jl:911
[5] _moi_get_result(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Gurobi.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, ::MathOptInterface.ConstraintDual, ::Vararg{Any})
@ JuMP C:\Users\Admin\.julia\packages\JuMP\ptoff\src\optimizer_interface.jl:640
[6] get(model::Model, attr::MathOptInterface.ConstraintDual, cr::ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape})
@ JuMP C:\Users\Admin\.julia\packages\JuMP\ptoff\src\optimizer_interface.jl:700
[7] _constraint_dual(con_ref::ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}, result::Int64)
@ JuMP C:\Users\Admin\.julia\packages\JuMP\ptoff\src\constraints.jl:1082
[8] #dual#20
@ C:\Users\Admin\.julia\packages\JuMP\ptoff\src\constraints.jl:1065 [inlined]
[9] dual
@ C:\Users\Admin\.julia\packages\JuMP\ptoff\src\constraints.jl:1061 [inlined]
[10] _broadcast_getindex_evalf
@ .\broadcast.jl:683 [inlined]
[11] _broadcast_getindex
@ .\broadcast.jl:656 [inlined]
[12] getindex
@ .\broadcast.jl:610 [inlined]
[13] macro expansion
@ .\broadcast.jl:974 [inlined]
[14] macro expansion
@ .\simdloop.jl:77 [inlined]
[15] copyto!
@ .\broadcast.jl:973 [inlined]
[16] copyto!
@ .\broadcast.jl:926 [inlined]
[17] copy
@ .\broadcast.jl:898 [inlined]
[18] materialize
@ .\broadcast.jl:873 [inlined]
[19] broadcast(f::typeof(dual), As::Vector{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}})
@ Base.Broadcast .\broadcast.jl:811
[20] broadcasted(#unused#::Base.Broadcast.ArrayStyle{JuMP.Containers.DenseAxisArray}, f::Function, args::JuMP.Containers.DenseAxisArray{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}, 1, Tuple{Vector{Int64}}, Tuple{JuMP.Containers._AxisLookup{Dict{Int64, Int64}}}})
@ JuMP.Containers C:\Users\Admin\.julia\packages\JuMP\ptoff\src\Containers\DenseAxisArray.jl:554
[21] broadcasted(::Function, ::JuMP.Containers.DenseAxisArray{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}, 1, Tuple{Vector{Int64}}, Tuple{JuMP.Containers._AxisLookup{Dict{Int64, Int64}}}})
@ Base.Broadcast .\broadcast.jl:1311
[22] DUC_Clearing(UnitCommitment_dt::UC_Data)
@ Main d:\Windows PC\IITB\Sem5\Optimization with Julia\Unit Commitment\Deterministic\DUC_3G.jl:130
[23] top-level scope
@ .\timing.jl:273