How can we find out the objective function value for all the points in the feasible region for a complex optimization problem using JuMP?
Is there a built-in method for doing so? If yes, then please specify the same.
If not, then please tell a way to do this.
There might be infinitely many points in the feasible region, how do you intend to “figure out the objection-function value” for such a large number of points?
I want to generate some random points let’s say 1000 or so.
I want to get a picture of how the objective function varies. Is there any extreme value?
The objective function now includes a weighted sum of the expected cost and the CVaR. For the input data set that I am using, I found that for only expected cost minimization the CVaR comes out to be very large, while is around the expected value for all other values between 0 and 1. So, I want to know is there a point in the feasible region which is resulting in such an extreme result.
For a problem, the feasible solutions can be infinite, If you want to get the extreme objective function values, you can just set the optimization aims to minimize the objective function or maximize it to get the corresponding solutions.
I want to visualize the objective function value for several random feasible points.
I assume that this just means your variable which represents CVaR is not “tight” and that it can take any value. This can happen because it is not in the objective function, and so there is no incentive for the solver to pick a small value.
If you have a reproducible example, we might be able to say more.
I want to visualize the objective function value for several random feasible points.
You can evaluate the objective function at a point as follows:
obj_f = objective_function(model)
point = Dict(
x => rand() # replace here with a feasible primal solution
for x in all_variables(model)
)
value(x -> point[x], obj_f)
It’s a lot harder to find a set of random feasible points.
This will also only be useful if your problem has two dimensions. What do you hope to visualize?
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")
When beta=1 the CVaR term does not appear in the objective, and therefore CVaR_Total_Cost can take any value. It is only guaranteed to be equal to CVaR of the distribution if it attains a minimum.
Just use a value of beta just less than 1, or manually compute the CVaR based on the solution of the costs.
Yes, that is true. CVaR can take any value for beta=1.
The value of CVaR for beta =1 is also obtained from the optimal solution. How do I manually compute the CVaR based on the solution of the costs?
If I want to limit its value then I have to add an upper bound on CVaR.
Perhaps to clarify in case it wasn’t clear, its not that the true CVaR of the distribution can take any value (although it can, provided the rest of the objective stays minimized), but that ζ and η can take any values, and so your expression for CVaR will not equal what you would get if you solved the problem without CVaR and then computed it based on the solution values.
It just doesn’t make sense to use this formulation of CVaR if it does not appear in the objective.