```
ESS = Model(Gurobi.Optimizer)
# Variables
Pgen = @variable(ESS, [bus_key, scenar], lower_bound=0, base_name ="Pgen")
PD = @variable(ESS, [scenar], lower_bound=0, base_name="PD")
PC = @variable(ESS, [scenar], lower_bound=0, base_name="PC")
Eb = @variable(ESS, [scenar], lower_bound=0, base_name="Eb")
max_Pb = @variable(ESS, lower_bound=0, base_name="max_Pb")
max_Eb = @variable(ESS, lower_bound=0, base_name="max_Eb")
PgenT = @variable(ESS, [scenar, bus_key], lower_bound=0, base_name ="PgenT")
#cost_genvar = @variable(ESS)
## constraints
Pgenlim = @constraint(ESS, [n=bus_key,s=sc], Pgen[n,s]<=PG[n,s], base_name="Pgenlim")
kcl = @constraint(ESS, [s=sc], sum(Pgen[n,s] for n in bus_key) + PD[s] - PC[s] - PL[s] == 0, base_name="kcl")
# Limits on Charging and Discharging Power, Energy, and Pmax, Emax
PDmax_up_lim =@constraint(ESS, [s=sc], PD[s] <= ESS_stat*max_Pb, base_name="PDmax_up_lim")
PCmax_up_lim =@constraint(ESS, [s=sc], PC[s] <= ESS_stat*max_Pb, base_name="PCmax_up_lim")
# Limits on Emax
Eb_up_lim = @constraint(ESS, [s=sc], Eb[s] <= ESS_stat*0.8*max_Eb, base_name="Eb_up_lim")
Eb_low_lim = @constraint(ESS, [s=sc], Eb[s] >= ESS_stat*0.2*max_Eb, base_name="Eb_low_lim")
Pbmax_up_lim = @constraint(ESS, max_Pb <= 10^6, base_name="Pbmax_up_lim")
# Limit on Discharge Duration
Ebmax_up_lim = @constraint(ESS, max_Eb <= DSmax*max_Pb/eta, base_name="Ebmax_up_lim")
Ebmax_low_lim = @constraint(ESS, max_Eb >= DSmin*max_Pb/eta, base_name="Ebmax_low_lim")
# Energy Balance constraint
Energy_balanceall = @constraint(ESS, [s=scenar_except_end], Eb[s+1] - Eb[s] + (PD[s]/eta - PC[s]*eta)*dur[s] == 0, base_name="Energy_balanceall")
Energy_balanceend = @constraint(ESS, [n_scen], Eb[1] - Eb[n_scen] + (PD[n_scen]/eta - PC[n_scen]*eta)*dur[n_scen] == 0, base_name="Energy_balanceend")
println(Energy_balanceend)
##
CRF=0.1130
for i in 1:n_bus
for j in 1:n_scen
@constraint(ESS, Pgen[i,j]==PgenT[j,i])
end
end
# cost_ESS
ESS_cost = CRF*(CapES*max_Eb*100 + CapPS*max_Pb*100) + FOM*max_Pb*100 + VOM*sum(PD[s]*dur[s] for s in scenar)*100*365
cost_genvar = sum(PgenT[j,i]*cov[i,k]*Pgen[k,j] for i in bus_key for j in scenar for k in bus_key)
if p == 0
@objective(ESS, Min, ESS_cost)
else
@objective(ESS, Min, ESS_cost + p*cost_genvar)
end
# solve
set_optimizer_attribute(ESS, "OutputFlag", 1)
set_optimizer_attribute(model, "FeasibilityTol", 1e-8)
println(cost_genvar)
```