using JuMP, Gurobi, Plots
function CVaR_3bus(CUD, M, α, β)
# Sets
w = ["w1", "w2", "w3", "w4", "w5", "w6"] # scenarios
d = ["d1"] # representative
t = ["t1", "t2", "t3"] # time periods in each representative
k = ["k1", "k2"] # groups for EV
g = ["g1", "g2", "g3"] # generators
s = ["s1", "s2", "s3"] # storage
b = ["b1", "b2", "b3"] # buses
l = ["l1", "l2", "l3", "l4", "l5", "l6"] # lines
T = collect(1:length(t)) # time periods in integer
#Lᴱ = filter(x -> parse(Int, x[2:end]) ≤ 3, l) # gives the existing lines [l1, l2, l3]
#Lᶜ = filter(x -> parse(Int, x[2:end]) > 3, l) # gives the candidate lines [l4, l5, l6]
Lᴱ = Dict("l1" => true,"l2" => true, "l3" => true,"l4" => false,"l5" => false, "l6" => false )
Lᶜ = Dict("l1" => false,"l2" => false, "l3" => false,"l4" => true,"l5" => true, "l6" => true )
Gᴰ = Dict("g1" => true, "g2" => true, "g3" => false) # Set of Dispatchable Units
Gᴵ = Dict("g1" => false, "g2" => false, "g3" => true) # Set of Intermittent Units
Gᵇ = Dict(
"b1"=>["g1"],
"b2"=>["g2"],
"b3"=>["g3"])
Sᵇ = Dict(
"b1"=>["s1"],
"b2"=>["s2"],
"b3"=>["s3"])
#=***********************************************************************************
Origin and Desitination Bus for Each Line
Lines Originating from bus and lines meeting to a destination bus
************************************************************************************=#
OL = Dict(
"l1" => "b1",
"l2" => "b2",
"l3" => "b3",
"l4" => "b1",
"l5" => "b2",
"l6" => "b3"
)
# Destination of Line 1
FL = Dict(
"l1" => "b2",
"l2" => "b3",
"l3" => "b1",
"l4" => "b2",
"l5" => "b3",
"l6" => "b1"
)
# Initialize an empty dictionary to store the result
L_OB = Dict{String, Vector{String}}()
# Iterate through the OL dictionary
for (lo, bo) in OL
# If the b key doesn't exist in the L_OB, create a new entry with an empty vector
if !haskey(L_OB, bo)
L_OB[bo] = []
end
# Add the l key to the corresponding b entry
push!(L_OB[bo], lo)
end
# Initialize an empty dictionary to store the result
L_DB = Dict{String, Vector{String}}()
# Iterate through the OL dictionary
for (ld, bd) in FL
# If the b key doesn't exist in the L_DB, create a new entry with an empty vector
if !haskey(L_DB, bd)
L_DB[bd] = []
end
# Add the l key to the corresponding b entry
push!(L_DB[bd], ld)
end
#=*****************************************************************************
Scalars
*****************************************************************************=#
#ηEV = 0.98 # Chargin and discharging efficiency of EV
#Eₑᴮm = 0.04 # Battery Capacity of EV
#γₑm = 0.05 # Minimum energy stored in EV
θᴰm = 1.5707963 # Maximum voltage angle deviation
#=*****************************************************************************
Parameters
*****************************************************************************=#
WD = Dict("d1" => 365*8) # Weight for each time period
γˢm = Dict("s1" => 0.1, "s2" => 0.1, "s3" => 0.1) # per unit energy that must remain in the storage
γₛ⁰ = Dict(
("s1","d1") => 0.7,
("s2","d1") => 0.7,
("s3","d1") => 0.7,) # Per unit energy in the first period(storage)
γₛEP = Dict("s1" => 6, "s2" => 6, "s3" => 6) # 6 hr discharge duration
ηS = Dict("s1" => 0.98, "s2" => 0.98, "s3" => 0.98) # Storage efficiency
# Define the table for Number of Electric Vehicles
EVₙ = Dict(
("k1", "b1") => Dict("w1" => 550, "w2" => 550, "w3" => 750, "w4" => 450, "w5" => 950, "w6" => 950),
("k1", "b2") => Dict("w1" => 800, "w2" => 800, "w3" => 1000, "w4" => 1000, "w5" => 1200, "w6" => 1200),
("k1", "b3") => Dict("w1" => 400, "w2" => 400, "w3" => 600, "w4" => 600, "w5" => 800, "w6" => 800),
("k2", "b1") => Dict("w1" => 300, "w2" => 300, "w3" => 500, "w4" => 500, "w5" => 700, "w6" => 750),
("k2", "b2") => Dict("w1" => 400, "w2" => 400, "w3" => 600, "w4" => 600, "w5" => 800, "w6" => 850),
("k2", "b3") => Dict("w1" => 150, "w2" => 150, "w3" => 350, "w4" => 350, "w5" => 550, "w6" => 550)
)
pi = Dict("w1" => 1/length(w), "w2" => 1/length(w), "w3" => 1/length(w), "w4" => 1/length(w), "w5" => 1/length(w), "w6" => 1/length(w))
Pᵢᴳm = Dict("g1" => 400, "g2" => 400, "g3" => 600) # Maximum Capacity of generting units
Eₛm = Dict("s1" => 100, "s2" => 100, "s3" => 100) # Maximum possible Energy capacity of the candidate storage
Pₛm = Dict("s1" => 100, "s2" => 100, "s3" => 100) # Maximum possible Power capacity of the candidate storage
PGᵁ = Dict("g1" => 0.25, "g2" => 0.25) # Per unit ramp up limit for the generating unit
PGᴰ = Dict("g1" => 0.25, "g2" => 0.25) # Per unit ramp down limit for the generating unit
Pₗm = Dict("l1" => 70, "l2" => 70, "l3" => 70, "l4" => 70, "l5" => 70, "l6" => 70) # Line Capacity
X = Dict("l1" => 0.01, "l2" => 0.01, "l3" => 0.01, "l4" => 0.01, "l5" => 0.01, "l6" => 0.01) # Line reactance
CIL = Dict("l4" => 200000, "l5" => 200000, "l6" => 200000) # Transmission Line Cost
#
# Table AD availability of intermittnet power units
AD = Dict(
("g3", "d1") => Dict("t1" => 0.4, "t2" => 0.9, "t3" => 0.6)
)
# Table Generation Investment Plan for Renewables
Pᵢᴳ = Dict(
("g1") => Dict("w1" => 300, "w2" => 150, "w3" => 300, "w4" => 150, "w5" => 300, "w6" => 150),
("g2") => Dict("w1" => 90, "w2" => 90, "w3" => 90, "w4" => 90, "w5" => 90, "w6" => 90),
("g3") => Dict("w1" => 100, "w2" => 300, "w3" => 100, "w4" => 300, "w5" => 100, "w6" => 300),
)
# Table Demand (MWh)
Pᴰ = Dict(
("b1","d1","t1") => Dict("w1" => 80, "w2" => 80, "w3" => 80, "w4" => 80, "w5" => 80, "w6" => 80),
("b1","d1","t2") => Dict("w1" => 100, "w2" => 100, "w3" => 100,"w4" => 100, "w5" => 100, "w6" => 100),
("b1","d1","t3") => Dict("w1" => 90, "w2" => 90, "w3" => 90, "w4" => 90, "w5" => 90, "w6" => 90),
("b2","d1","t1") => Dict("w1" => 220, "w2" => 220, "w3" => 220,"w4" => 220, "w5" => 220, "w6" => 220),
("b2","d1","t2") => Dict("w1" => 250, "w2" => 250, "w3" => 250,"w4" => 250, "w5" => 250, "w6" => 250),
("b2","d1","t3") => Dict("w1" => 230, "w2" => 230, "w3" => 230,"w4" => 230, "w5" => 230, "w6" => 230),
("b3","d1","t1") => Dict("w1" => 35, "w2" => 35, "w3" => 35, "w4" => 35, "w5" => 35, "w6" => 35),
("b3","d1","t2") => Dict("w1" => 50, "w2" => 50, "w3" => 50, "w4" => 50, "w5" => 50, "w6" => 50),
("b3","d1","t3") => Dict("w1" => 40, "w2" => 40, "w3" => 40, "w4" => 40, "w5" => 40, "w6" => 40)
)
# Table PD
PEV_pu = Dict(
("k1","b1","d1","t1") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0),
("k1","b1","d1","t2") => Dict("w1" => 15.00, "w2" => 15.00, "w3" => 15.00, "w4" => 15.00, "w5" => 15.00, "w6" => 15.00),
("k1","b1","d1","t3") => Dict("w1" => 5.408, "w2" => 5.408, "w3" => 5.408, "w4" => 5.408, "w5" => 5.408, "w6" => 5.408),
("k1","b2","d1","t1") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0),
("k1","b2","d1","t2") => Dict("w1" => 15.00, "w2" => 15.00, "w3" => 15.00, "w4" => 15.00, "w5" => 15.00, "w6" => 15.00),
("k1","b2","d1","t3") => Dict("w1" => 5.408, "w2" => 5.408, "w3" => 5.408, "w4" => 5.408, "w5" => 5.408, "w6" => 5.408),
("k1","b3","d1","t1") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0),
("k1","b3","d1","t2") => Dict("w1" => 15.00, "w2" => 15.00, "w3" => 15.00, "w4" => 15.00, "w5" => 15.00, "w6" => 15.00),
("k1","b3","d1","t3") => Dict("w1" => 5.408, "w2" => 5.408, "w3" => 5.408, "w4" => 5.408, "w5" => 5.408, "w6" => 5.408),
("k2","b1","d1","t1") => Dict("w1" => 15.00, "w2" => 15.00, "w3" => 15.00, "w4" => 15.00, "w5" => 15.00, "w6" => 15.00),
("k2","b1","d1","t2") => Dict("w1" => 5.408, "w2" => 5.408, "w3" => 5.408, "w4" => 5.408, "w5" => 5.408, "w6" => 5.408),
("k2","b1","d1","t3") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0),
("k2","b2","d1","t1") => Dict("w1" => 15.00, "w2" => 15.00, "w3" => 15.00, "w4" => 15.00, "w5" => 15.00, "w6" => 15.00),
("k2","b2","d1","t2") => Dict("w1" => 5.408, "w2" => 5.408, "w3" => 5.408, "w4" => 5.408, "w5" => 5.408, "w6" => 5.408),
("k2","b2","d1","t3") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0),
("k2","b3","d1","t1") => Dict("w1" => 15.00, "w2" => 15.00, "w3" => 15.00, "w4" => 15.00, "w5" => 15.00, "w6" => 15.00),
("k2","b3","d1","t2") => Dict("w1" => 5.408, "w2" => 5.408, "w3" => 5.408, "w4" => 5.408, "w5" => 5.408, "w6" => 5.408),
("k2","b3","d1","t3") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0)
)
# Initialize a new dictionary to store the multiplied values
PEV = Dict{Tuple, Dict}()
# Iterate through the keys of EVₙ
for (EV_key, EV_value) in EVₙ
# Extract "k" and "b" values from the EV_key
kEV, bEV = EV_key
# Iterate through the keys of PD
for (PEV_pu_key, PEV_pu_value) in PEV_pu
# Check if the "k" and "b" values match
if (kEV, bEV) == (PEV_pu_key[1], PEV_pu_key[2])
# Multiply the corresponding values
multiplied_dict = Dict{String, Float64}()
for (wEV, EV_val) in EV_value
multiplied_dict[wEV] = EV_val * PEV_pu_value[wEV]/1000
end
# Store the result in the result_dict
PEV[PEV_pu_key] = multiplied_dict
end
end
end
# Table Annualized Energy Investment Cost
CEᵢˢ = Dict(
("s1") => Dict("w1" => 500, "w2" => 500, "w3" => 500, "w4" => 500, "w5" => 500, "w6" => 500),
("s2") => Dict("w1" => 500, "w2" => 500, "w3" => 500, "w4" => 500, "w5" => 500, "w6" => 500),
("s3") => Dict("w1" => 500, "w2" => 500, "w3" => 500, "w4" => 500, "w5" => 500, "w6" => 500),
)
# Table Annualized Power Investment Cost
CPᵢˢ = Dict(
("s1") => Dict("w1" => 3000, "w2" => 3000, "w3" => 3000, "w4" =>3000, "w5" =>3000, "w6" =>3000),
("s2") => Dict("w1" => 3000, "w2" => 3000, "w3" => 3000, "w4" =>3000, "w5" =>3000, "w6" =>3000),
("s3") => Dict("w1" => 3000, "w2" => 3000, "w3" => 3000, "w4" =>3000, "w5" =>3000, "w6" =>3000)
)
# Table 'Operating Cost (generating Units)'
CᴳD = Dict(
("g1") => Dict("w1" => 0.250, "w2" => 0.250, "w3" => 0.250, "w4" => 0.250, "w5" => 0.250, "w6" => 0.250),
("g2") => Dict("w1" => 0.400, "w2" => 0.400, "w3" => 0.400, "w4" => 0.400, "w5" => 0.400, "w6" => 0.400),
("g3") => Dict("w1" => 0.000, "w2" => 0.000, "w3" => 0.000, "w4" => 0.000, "w5" => 0.000, "w6" => 0.000)
)
# Table Discharge Cost(Storage)
CˢD = Dict(
("s1") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0),
("s2") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0),
("s3") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0)
)
# Table CSC
CˢC = Dict(
("s1") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0),
("s2") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0),
("s3") => Dict("w1" => 0, "w2" => 0, "w3" => 0, "w4" => 0, "w5" => 0, "w6" => 0)
)
#=***********************************************************************************
Decalaration of Variables and Mathematical Characterization
*************************************************************************************=#
# Create a JuMP model with Gurobi solver
model = Model(Gurobi.Optimizer)
set_attribute(model, "OutputFlag", 1)
set_attribute(model, "FeasibilityTol", 1e-8)
set_attribute(model, "MIPGap", 1e-8)
set_attribute(model, "MIPGapAbs", 0)
# Variables
@variables(model, begin
# First stage variables
eᵢˢ[s in s] ≥ 0 # Storage energy capacity installed
pᵢˢ[s in s] ≥ 0 # Storage power capacity installed
yᴸ[l[4:end]], Bin # Binary variable that is equal to 1 if line l is built
# Second Stage Variable
pᴳD[g in g, d in d, t in t, w in w] ≥ 0
pᴳDS[g in g, d in d, t in t, w in w] ≥ 0 # Intermittent power spillage
pˢC[s in s, d in d, t in t, w in w] ≥ 0 # Power charged from storages
pˢD[s in s, d in d, t in t, w in w] ≥ 0 # Power discharged from storages
eˢ[s in s, d in d, t in t, w in w] ≥ 0 # Energy stored in storages
pᴸD[l in l, d in d, t in t, w in w] # Power flows
θᴰ[b in b, d in d, t in t, w in w] # Voltage angles
pᵁD[b in b, d in d, t in t, w in w] ≥ 0 # Unserved demand
# CVaR Formulation Variables
η[w in w] ≥0
ζ
end)
#=***********************************************************************************
Constraints and Equations
*************************************************************************************=#
@expressions(model,begin
# Expressions for the Energy Balance Cosntraints
Gen[b in b, d in d, t in t, w in w], sum(pᴳD[gb, d, t, w] for gb in Gᵇ[b]) # Sum of all generation at a bus
ESSgen[b in b, d in d, t in t, w in w], sum(pˢD[sb, d, t, w] - pˢC[sb, d, t, w] for sb in Sᵇ[b]) # Sum of difference of discharge and cahrge power for all storages at a bus
PF_final_bus[b in b, d in d, t in t, w in w], sum(pᴸD[lfb, d, t, w] for lfb in L_DB[b]) # sum of power flows for lines originating from each bus
PF_origin_bus[b in b, d in d, t in t, w in w], sum(pᴸD[lob, d, t, w] for lob in L_OB[b]) # sum of power flows for lines meeting to each bus
Demand[b in b, d in d, t in t, w in w], Pᴰ[b, d, t][w] # Demand at each bus
EVdemand[b in b, d in d, t in t, w in w], sum(PEV[k, b, d, t][w] for k in k) # sum of EV demand for each group at each bus
DUnserved[b in b, d in d, t in t, w in w], pᵁD[b, d, t, w] # load curtailed for each bus, dat, time, scenario
# Expressions for the Objective Function
#Cost_Lᶜ, sum(pi[w]*sum([yᴸ[l[i]] * CIL[l[i]] for i in collect(1:length(l)) if i>3]) for w in w) # Annualized Cost of Candidate Transmission Lines
Cost_Lᶜ, sum([yᴸ[l[i]] * CIL[l[i]] for i in collect(1:length(l)) if i>3])
CapCost_ESS, sum(pi[w]*sum([CEᵢˢ[es][w] * eᵢˢ[es]+ CPᵢˢ[es][w]*pᵢˢ[es] for es in s]) for w in w) # Annualized Energy and Power cost of Energ Storage
OpCost_Gen, sum(pi[w]*sum(sum(WD[d]*sum(CᴳD[gd][w]*pᴳD[gd, d, t, w] for gd in g) for t in t) for d in d) for w in w)
OpCost_ESS, sum(pi[w]*sum(sum(WD[d]*sum(CˢD[sb][w]*pˢD[sb, d, t, w] + CˢC[sb][w]*pˢC[sb, d, t, w] for sb in s) for t in t) for d in d) for w in w)
Cost_DUnserved, sum(pi[w]*sum(sum(WD[d]*sum(CUD*pᵁD[b, d, t, w] for b in b) for t in t) for d in d) for w in w)
#=
OpCost_Gen, sum(pi[w]*sum(sum(WD[d]*sum(CᴳD[gd]["w1"]*pᴳD[gd, d, t, w] for gd in g)+
sum(CˢD[sb][w]*pˢD[sb, d, t, w] + CˢC[sb][w]*pˢC[sb, d, t, w] for sb in s)+
sum(CUD*pᵁD[b, d, t, w] for b in b) for t in t) for d in d) for w in w)
OpCost_Gen, sum(pi[w]*sum(sum(WD[d]*sum(CᴳD[gd]["w1"]*pᴳD[gd, d, t, w] for gd in g)+
sum(CˢD[sb][w]*pˢD[sb, d, t, w] + CˢC[sb][w]*pˢC[sb, d, t, w] for sb in s)+
sum(CUD*pᵁD[b, d, t, w] for b in b) for t in t) for d in d) for w in w)
=#
CVaR_Total_Cost, ζ + (1/(1-α))*(pi[w]*sum(η[w] for w in w))
end)
@constraints(model, begin
# Investment Constraints
InvestESS_01[s in s], eᵢˢ[s] ≤ Eₛm[s] # Energy capcaity ≤ max possible value
InvestESS_02[s in s], pᵢˢ[s] ≤ Pₛm[s] # Power capcaity ≤ max possible value
InvestESS_03[s in s], eᵢˢ[s] ≥ γₛEP[s] * pᵢˢ[s] # defining discharge duration for the BESS
# Generating Units Constraints
Gen_01[g in g, d in d, t in t, w in w; Gᴰ[g]], pᴳD[g, d, t, w] ≤ Pᵢᴳ[g][w] # Dispatchable gen ≤ Capacity
Gen_02[g in g, d in d, t in t, w in w; Gᴵ[g]], pᴳD[g, d, t, w] + pᴳDS[g, d, t, w] ≤ AD[g, d][t] * Pᵢᴳ[g][w] # RES Gen + RES Spillage ≤ AF*Capacity
Gen_03[g in g, d in d, i in T, w in w; Gᴰ[g] && i!=1], pᴳD[g, d, t[i], w] - pᴳD[g, d, t[i-1], w] ≤ PGᵁ[g] * Pᵢᴳ[g][w] # Ramp up Constraints
Gen_04[g in g, d in d, i in T, w in w; Gᴰ[g] && i!=1], pᴳD[g, d, t[i-1], w] - pᴳD[g, d, t[i], w] ≤ PGᴰ[g] * Pᵢᴳ[g][w] # Remap Down Constraints
# Storage Unit Constraints
ESS_01[s in s, d in d, t in t, w in w], pˢD[s, d, t, w] ≤ pᵢˢ[s] # discharging power ≤ Power capacity
ESS_02[s in s, d in d, t in t, w in w], pˢC[s, d, t, w] ≤ pᵢˢ[s] # charging power ≤ Power capacity
ESS_03[s in s, d in d, i in T, w in w; i!=1], eˢ[s, d, t[i], w] == eˢ[s, d, t[i-1], w] + ηS[s] * pˢC[s, d, t[i], w] - (1 / ηS[s]) * pˢD[s, d, t[i], w] # Energy balance for time period greater than 1
ESS_04[s in s, d in d, i in T, w in w; i==1], eˢ[s, d, t[i], w] == γₛ⁰[s, d] * eᵢˢ[s]+ ηS[s] * pˢC[s, d, t[i], w] - (1 / ηS[s]) * pˢD[s, d, t[i], w] # Energy balance for first time period
ESS_05[s in s, d in d, t in t, w in w], eˢ[s, d, t, w] >= γˢm[s] * eᵢˢ[s] # energy at any time ≥ min value possible
ESS_06[s in s, d in d, t in t, w in w], eˢ[s, d, t, w] <= eᵢˢ[s] # energy at any time ≤ max value possible
ESS_07[s in s, d in d, t in t, w in w; t == "t3"], eˢ[s, d, t, w] ≥ γₛ⁰[s, d] * eᵢˢ[s] # SOC at the end of the day must be greater than a specified value
# Transmission Constraints
TL_01[l in l, d in d, t in t, w in w; Lᴱ[l]], pᴸD[l, d, t, w] == 1 / X[l] * (θᴰ[OL[l], d, t, w] - θᴰ[FL[l], d, t, w] )
TL_02[l in l, d in d, t in t, w in w; Lᶜ[l]], pᴸD[l, d, t, w] ≥ 1 / X[l] * (θᴰ[OL[l], d, t, w] - θᴰ[FL[l], d, t, w] )-(1-yᴸ[l])*M
TL_03[l in l, d in d, t in t, w in w; Lᶜ[l]], pᴸD[l, d, t, w] ≤ 1 / X[l] * (θᴰ[OL[l], d, t, w] - θᴰ[FL[l], d, t, w] )+(1-yᴸ[l])*M
TL_04[l in l, d in d, t in t, w in w; Lᴱ[l]], pᴸD[l, d, t, w] ≤ Pₗm[l]
TL_05[l in l, d in d, t in t, w in w; Lᴱ[l]], pᴸD[l, d, t, w] ≥ -Pₗm[l]
TL_06[l in l, d in d, t in t, w in w; Lᶜ[l]], pᴸD[l, d, t, w] ≤ Pₗm[l] * yᴸ[l]
TL_07[l in l, d in d, t in t, w in w; Lᶜ[l]], pᴸD[l, d, t, w] ≥ -Pₗm[l] * yᴸ[l]
TL_08[b in b, d in d, t in t, w in w], θᴰ[b, d, t, w] ≥ -θᴰm
TL_09[b in b, d in d, t in t, w in w], θᴰ[b, d, t, w] ≤ θᴰm
# Energy Balance Constraints
EnergyBalance01[b in b, d in d, t in t, w in w],
Gen[b,d,t,w] + ESSgen[b,d,t,w] - PF_origin_bus[b,d,t,w] + PF_final_bus[b,d,t,w] == Demand[b,d,t,w] + EVdemand[b,d,t,w] - DUnserved[b,d,t,w]
EnergyBalance02[b in b, d in d, t in t, w in w], DUnserved[b,d,t,w] ≤ Demand[b,d,t,w] + EVdemand[b,d,t,w] # unserved demand ≤ the inflexible demand
# Constraints Related to CVaR for each scenario
CVaR[w in w], (sum([yᴸ[l[i]] * CIL[l[i]] for i in collect(1:length(l)) if i>3]) +
sum([CEᵢˢ[es][w] * eᵢˢ[es]+ CPᵢˢ[es][w]*pᵢˢ[es] for es in s])+
sum(sum(WD[d]*sum(CᴳD[gd][w]*pᴳD[gd, d, t, w] for gd in g) for t in t) for d in d)+
sum(sum(WD[d]*sum(CˢD[sb][w]*pˢD[sb, d, t, w] + CˢC[sb][w]*pˢC[sb, d, t, w] for sb in s) for t in t) for d in d)+
sum(sum(WD[d]*sum(CUD*pᵁD[b, d, t, w] for b in b) for t in t) for d in d)) -ζ ≤ η[w]
end)
# Objective function
#z= @objective(model, Min, 0.001*(Cost_Lᶜ + CapCost_ESS + OpCost_Gen)) # 0.0001 to convert € to k€
z= @objective(model, Min, 0.001*(β*(Cost_Lᶜ + CapCost_ESS + OpCost_Gen + OpCost_ESS + Cost_DUnserved)+(1-β)*CVaR_Total_Cost))
# Solve the model
optimize!(model)
#=****************************************************************************************************
Results
****************************************************************************************************=#
#=
println("Investment Cost ")
println("Storage Investment Cost: ", 0.001*value(CapCost_ESS), " k€")
println("Transmission Line Investment Cost: ", 0.001*value(Cost_Lᶜ), " k€")
println("Generation Operational Cost: ", 0.001*(value(OpCost_Gen)), " k€" )
println("Cost of Demand Unserved: ", 0.001*(value(Cost_DUnserved)), " k€" )
println("Candidate Transmission Lines:", value.(yᴸ))
println("Storage Requirement in MWh and MW ")
println([value.(eᵢˢ) value.(pᵢˢ)])
println(value.(pᴸD))
println(value.(pˢC))
println(value.(pˢD))
println(0.001*value(CVaR_Total_Cost))
println("Total Cost: ", (0.001*(value(Cost_Lᶜ + CapCost_ESS + OpCost_Gen + OpCost_ESS + Cost_DUnserved))), " k€" )
=#
Inv_Cost = value(Cost_Lᶜ + CapCost_ESS)*0.001
Ex_Op_Cost = value(OpCost_Gen )*0.001
Ex_DUnserved = value(Cost_DUnserved)*0.001
Ex_Tot_Cost = value(Cost_Lᶜ + CapCost_ESS + OpCost_Gen + OpCost_ESS + Cost_DUnserved).*0.001
CVaR_Cost = value(CVaR_Total_Cost)*0.001
#=
println("Invextment Cost(k€): ", )
println("Expected Operation Cost(k€): ", )
println("Expected Demand Unserved Cost(k€): ", )
println("Expected Total Cost: ", )
println("CVaR :", )
=#
return Inv_Cost, Ex_Op_Cost, Ex_DUnserved, Ex_Tot_Cost, CVaR_Cost
end
CUD = 1000 # Unreserves Demand Cost
M = 1000 # Big M
α = 0.90 # confidence Level for CVaR formulation
β = 0:0.01:1 # Weighting Factor for Expected Cost and Risk associated with it
Tot_Cost = Vector{Float64}(undef, length(β))
CVaR_Cost = Vector{Float64}(undef, length(β))
for i = 1:length(β)
IV, OC, DC, TC, CC = CVaR_3bus(CUD, M, α, β[i])
Tot_Cost[i] = TC
CVaR_Cost[i] = CC
end
scatter(CVaR_Cost, Tot_Cost, zcolor=β, color=:blues, xlabel="CVaR (k€)", ylabel="Expected Cost (k€)", legend=false, title="Efficient Frontier for α=0.9")