Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

FIx depreciation warnings + more #4

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

FIx depreciation warnings + more #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4e8b2d2
FIx depreciation warnings
Sep 21, 2016
7a6ceb0
Swapped direction of chance constraints
Oct 3, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@

*.mat

113 changes: 53 additions & 60 deletions codejl/ccopfmodel_simulation.jl
100644 → 100755
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,28 +30,28 @@ function solve_ccopf(logfilename, refbus, line_probability_threshold, gen_probab
numlines = length(lines) # Number of lines
numfarms = length(farms) # Number of farms

const Γ=iround(robust_budget*numfarms) # Budget of uncertainty
const Γ=round(Integer,robust_budget*numfarms) # Budget of uncertainty
const VoLL = 10_000 # Value of lost load

# Define the name of the model and solver settings
m = ChanceModel(solver=GurobiSolver(Method=1,BarHomogeneous=1))

# Define variables:
# Define variables:
# Non-negative hourly average power output bounded by the maximum limit on each generator:
@defVar(m, 0 <= pbar[i=generatorlist] <= buses[i].Pgmax)
@variable(m, 0 <= pbar[i=generatorlist] <= buses[i].Pgmax)
# Regulation participation factor between 0 and 1:
@defVar(m, 0 <= alpha[i=generatorlist] <= ((buses[i].Pgmax == 0) ? 0 : 1))
@variable(m, 0 <= alpha[i=generatorlist] <= ((buses[i].Pgmax == 0) ? 0 : 1))
# Power flow in each transmission line should be within the limits:
@defVar(m, -lines[i].u <= barf[i=1:numlines] <= lines[i].u)
@variable(m, -lines[i].u <= barf[i=1:numlines] <= lines[i].u)

@defVar(m, delta[1:(numbuses-1)])
@defVar(m, h[1:(numbuses-1)]) # \hat B^-1 \bar p
@defVar(m, θ[1:(numbuses-1)]) # thetahat
@variable(m, delta[1:(numbuses-1)])
@variable(m, h[1:(numbuses-1)]) # \hat B^-1 \bar p
@variable(m, θ[1:(numbuses-1)]) # thetahat

@defVar(m, slack_bus[1:(numbuses-1)] >= 0)
@variable(m, slack_bus[1:(numbuses-1)] >= 0)

# Define random variable for each wind farm with a given mean and variance
@defIndepNormal(m, winderr[i=1:numfarms], mean=farms[i].deviation_mean_interval, var=farms[i].deviation_sigmasq_interval)
@indepnormal(m, winderr[i=1:numfarms], mean=farms[i].deviation_mean_interval, var=farms[i].deviation_sigmasq_interval)

# Set up the reference bus:
@assert refbus == numbuses
Expand All @@ -60,42 +60,42 @@ function solve_ccopf(logfilename, refbus, line_probability_threshold, gen_probab
# Define constraints on participation factors:
for i in 2:length(hydro_idx)
# Enforce uniform participation factors on hydro power plants
@addConstraint(m, alpha[hydro_idx[i-1]] == alpha[hydro_idx[i]])
@constraint(m, alpha[hydro_idx[i-1]] == alpha[hydro_idx[i]])
end
# Set participation factors of nuclear generators to 0 (provide no regulation due to the ramping/cycling limitations)
@addConstraint(m, nuclear_alpha_fixed[i=nuclear_idx], alpha[i] == 0)
@constraint(m, nuclear_alpha_fixed[i=nuclear_idx], alpha[i] == 0)
# Integrality constraint on participation factors:
@addConstraint(m, sumalpha, sum(alpha) == 1)
@constraint(m, sumalpha, sum(alpha) == 1)

# Exogenous outputs of hydro and nuclear generators:
# Enforce the hourly average power output (exogenous value, obtained from a separate hydro scheduling procedure) of wind farms
@addConstraint(m, hydro_ub[k=1:length(hydro_idx)], pbar[hydro_idx[k]] <= hydro_limit[k])
# Enforce the hourly average power output of nuclear power plants (always committed and output at their maximum power output - "must run" units)
@addConstraint(m, nuclear_fixed[i=nuclear_idx], pbar[i] == buses[i].Pgmax)
# Enforce the hourly average power output (exogenous value, obtained from a separate hydro scheduling procedure) of wind farms
@constraint(m, hydro_ub[k=1:length(hydro_idx)], pbar[hydro_idx[k]] <= hydro_limit[k])
# Enforce the hourly average power output of nuclear power plants (always committed and output at their maximum power output - "must run" units)
@constraint(m, nuclear_fixed[i=nuclear_idx], pbar[i] == buses[i].Pgmax)

# System-wide power balance constraint:
# System-wide power balance constraint:
sumload = sum([getLoad(b) for b in buses]) # system-wide load
sumloadminusmeanwind = sum([getLoadMinusMeanWind(b) for b in buses]) # system-wide net load, i.e. the system-wide load minus the system-wide wind
# Constraint on power balance:
@addConstraint(m, balance, sum(pbar) == sumloadminusmeanwind)
@constraint(m, balance, sum(pbar) == sumloadminusmeanwind)

Brow = Bhatsp'

# \hat B delta = \alpha
@addConstraint(m, defalpha[i=1:(numbuses-1)],
@constraint(m, defalpha[i=1:(numbuses-1)],
sum{Brow.nzval[idx]*delta[Brow.rowval[idx]], idx in Brow.colptr[i]:(Brow.colptr[i+1]-1)} -
(isgen(buses[i]) ? alpha[i] : 0) == 0)

# \hat B h = p
@addConstraint(m, defh[i=1:(numbuses-1)],
@constraint(m, defh[i=1:(numbuses-1)],
sum{Brow.nzval[idx]*h[Brow.rowval[idx]], idx in Brow.colptr[i]:(Brow.colptr[i+1]-1)} -
(isgen(buses[i]) ? pbar[i] : 0) - slack_bus[i] == 0)

# theta = \hat B^{-1}(w - d) + h = \hat B^{-1}(p + w - d)
@addConstraint(m, BB[i=1:(numbuses-1)], θ[i]-h[i] == Binvb[i])
@constraint(m, BB[i=1:(numbuses-1)], θ[i]-h[i] == Binvb[i])

# Calculating flows in transmission lines:
@addConstraint(m, flowangle[i=1:numlines], barf[i] +
@constraint(m, flowangle[i=1:numlines], barf[i] +
(lines[i].tail != numbuses ? -lines[i].y*θ[lines[i].tail] : 0) +
(lines[i].head != numbuses ? lines[i].y*θ[lines[i].head] : 0) == 0)

Expand All @@ -117,39 +117,38 @@ function solve_ccopf(logfilename, refbus, line_probability_threshold, gen_probab
ccexpr += fexpr*winderr[j]
end
prob_threshold = (i in loaded_lines_idx) ? loaded_lines_probability_threshold : line_probability_threshold

# Formulate chance constraints for the + and - transmission limits:
addConstraint(m, barf[i] + lines[i].y*ccexpr >= getThermalCapacity(lines[i], mvaBase)*(1-linebufferamt), with_probability=prob_threshold, uncertainty_budget_mean=Γ, uncertainty_budget_variance=Γ)
addConstraint(m, barf[i] + lines[i].y*ccexpr <= -getThermalCapacity(lines[i], mvaBase), with_probability= prob_threshold, uncertainty_budget_mean=Γ, uncertainty_budget_variance=Γ)
@constraint(m, barf[i] + lines[i].y*ccexpr <= getThermalCapacity(lines[i], mvaBase)*(1-linebufferamt), with_probability=1-prob_threshold, uncertainty_budget_mean=Γ, uncertainty_budget_variance=Γ)
@constraint(m, barf[i] + lines[i].y*ccexpr >= -getThermalCapacity(lines[i], mvaBase), with_probability= 1-prob_threshold, uncertainty_budget_mean=Γ, uncertainty_budget_variance=Γ)
end

# Chance constraints for the maximum output of generators:
sum_deviations = sum([winderr[i] for i in 1:numfarms])
for i in generatorlist
addConstraint(m, pbar[i] - sum_deviations*alpha[i] >= buses[i].Pgmax, with_probability=gen_probability_threshold, uncertainty_budget_mean=Γ, uncertainty_budget_variance=Γ)
@constraint(m, pbar[i] - sum_deviations*alpha[i] <= buses[i].Pgmax, with_probability=1-gen_probability_threshold, uncertainty_budget_mean=Γ, uncertainty_budget_variance=Γ)
end

# Chance constraints on the upward (+) and downward(-) ramping:
for (bus_idx, rampmax) in ramp
addConstraint(m, -sum_deviations*alpha[bus_idx] >= rampmax, with_probability=gen_probability_threshold, uncertainty_budget_mean=Γ, uncertainty_budget_variance=Γ)
addConstraint(m, -sum_deviations*alpha[bus_idx] <= -rampmax, with_probability=gen_probability_threshold, uncertainty_budget_mean=Γ, uncertainty_budget_variance=Γ)
@constraint(m, -sum_deviations*alpha[bus_idx] <= rampmax, with_probability=1-gen_probability_threshold, uncertainty_budget_mean=Γ, uncertainty_budget_variance=Γ)
@constraint(m, -sum_deviations*alpha[bus_idx] >= -rampmax, with_probability=1-gen_probability_threshold, uncertainty_budget_mean=Γ, uncertainty_budget_variance=Γ)
end


# Set the objective function
@setObjective(m, Min, sum{buses[i].pi1*pbar[i]+buses[i].pi2*(pbar[i]^2 + sumvar*alpha[i]^2), i in generatorlist} + VoLL*sum(slack_bus))
# Set the objective function
@objective(m, Min, sum{buses[i].pi1*pbar[i]+buses[i].pi2*(pbar[i]^2 + sumvar*alpha[i]^2), i in generatorlist} + VoLL*sum(slack_bus))



tic() # track the execution time
status = solvechance(m,method=:Cuts, debug=false) # Solve the model
status = solve(m,method=:Cuts, debug=false) # Solve the model
solvetime = toq() # record the execution time
if status != :Optimal
# if the model isn't solved optimally, return NaNs
return status, NaN, solvetime, fill(NaN, numbuses), fill(NaN, numbuses), fill(NaN, numlines), fill(NaN, numbuses-1)
end
# Return the optimal values if the model is solved optimally
return status, getObjectiveValue(m), solvetime, getValue(alpha), getValue(pbar), getValue(barf), getValue(slack_bus)
return status, getobjectivevalue(m), solvetime, getvalue(alpha), getvalue(pbar), getvalue(barf), getvalue(slack_bus)

end

Expand All @@ -160,11 +159,11 @@ logfilename, casefilename, windfilename, costsfilename, refbus, line_probability
numbuses, numgens, generatorlist, numbranches, buses, lines = readcase(casefilename, costsfilename, refbus, loadscale, thermalLimitscale)

# Initialize index of specific (hydro, nucs, coal, gas) generators and transmission lines (loaded)
hydro_idx = int(readcsv(extras["hydro_indices"])[:])
nuclear_idx = int(readcsv(extras["nuclear_indices"])[:])
coal_idx = int(readcsv(extras["coal_indices"])[:])
gas_idx = int(readcsv(extras["gas_indices"])[:])
loaded_lines_idx = int(readcsv(extras["loaded_lines_file"])[:])
hydro_idx = round(Int64,readcsv(extras["hydro_indices"])[:])
nuclear_idx = round(Int64,readcsv(extras["nuclear_indices"])[:])
coal_idx = round(Int64,readcsv(extras["coal_indices"])[:])
gas_idx = round(Int64,readcsv(extras["gas_indices"])[:])
loaded_lines_idx = round(Int64,readcsv(extras["loaded_lines_file"])[:])
# Threshold for the loaded line:
loaded_lines_probability_threshold = float(extras["loaded_lines_probability_threshold"])

Expand All @@ -175,7 +174,7 @@ loaded_lines_probability_threshold = float(extras["loaded_lines_probability_thre
@assert issubset(gas_idx, generatorlist)

# Define the set of indices for hydro generators. Will be used to enforce the uniform participation factors.
fixed_alpha_subsets = [hydro_idx]
fixed_alpha_subsets = collect(hydro_idx)

# Extract the names of additional input files from the outputs of redconfig.jl (line 142):
wind_forecast_file = extras["wind_forecast_file"]
Expand All @@ -190,7 +189,7 @@ wind_sigma = readcsv(wind_sigma_file)
loads = readcsv(load_file)
hydro_limit = readcsv(hydro_limit_file)
ramp_csv = readcsv(ramp_file)
ramp = [(int(ramp_csv[i,1]),ramp_csv[i,2]) for i in 1:size(ramp_csv,1)]
ramp = [(round(Int,ramp_csv[i,1]),ramp_csv[i,2]) for i in 1:size(ramp_csv,1)]

# Extract the name of the out file from the outputs of redconfig.jl (line 142):
output_matfile = extras["output_matfile"]
Expand All @@ -200,7 +199,7 @@ robust_mean_lower_file = extras["robust_mean_lower_file"]
robust_mean_upper_file = extras["robust_mean_upper_file"]
robust_sigma_lower_file = extras["robust_sigma_lower_file"]
robust_sigma_upper_file = extras["robust_sigma_upper_file"]
# Read the input range on the means and sigma
# Read the input range on the means and sigma
robust_mean_lower = readcsv(robust_mean_lower_file)
robust_mean_upper = readcsv(robust_mean_upper_file)
robust_sigma_lower = readcsv(robust_sigma_lower_file)
Expand Down Expand Up @@ -236,7 +235,7 @@ numfarms = size(wind_forecast,1)

farms = Farm[]
for i in 1:numfarms
busidx = int(wind_forecast[i,1])
busidx = round(Int,wind_forecast[i,1])
f = Farm(busidx, 0, 0)
f.colbarinv = Array(Float64, numbuses)
push!(farms, f)
Expand All @@ -261,23 +260,23 @@ end
# Arrays for storing solution values
objvals = Float64[]
solvetimes = Float64[]
statuses = ASCIIString[]
alphavals = {}
pbarvals = {}
barfvals = {}
slackvals = {}
statuses = String[]
alphavals = []
pbarvals = []
barfvals = []
slackvals = []

for t in 1:numTimeSteps
println("\n\n######## T = $t\n")
for i in 1:size(loads,1)
busidx = int(loads[i,1])
busidx = round(Int,loads[i,1])
Pd = float(loads[i,t+1])
buses[busidx].Pd = Pd
end
for i in 1:numfarms
busidx = int(wind_forecast[i,1])
busidx = round(Int,wind_forecast[i,1])
mean = loadscale*float(wind_forecast[i,t+1])

std = loadscale*float(wind_sigma[i,t+1])
farms[i].mean = mean
farms[i].stddev = std
Expand All @@ -303,7 +302,7 @@ for t in 1:numTimeSteps
for i in 1:numbuses
setMeanWindMinusLoad(buses[i], farms)
end

sumvar = 0.0
for f in farms
sumvar += f.stddev^2
Expand All @@ -314,12 +313,12 @@ for t in 1:numTimeSteps
Binvb = zeros(numbuses)
Binvb[1:numbuses .!= refbus] = Blu\barb


# Solve the ccopf model
status, objval, solvetime, alphaval, pbarval, barfval, slackval =
solve_ccopf(logfilename, refbus, line_probability_threshold, gen_probability_threshold, linebufferamt, mvaBase, loadscale, thermalLimitscale, buses, lines, farms, generatorlist, Binvb, Bhatsp, sumvar, hydro_idx, nuclear_idx, hydro_limit[:,t+1], ramp, robust_budget, loaded_lines_idx, loaded_lines_probability_threshold)

# Write the optimal solution:

# Write the optimal solution:
push!(objvals, objval)
push!(solvetimes, solvetime)
push!(statuses, string(status))
Expand Down Expand Up @@ -364,9 +363,3 @@ write(mfile, "pbarvals", pbarmat)
write(mfile, "barfvals", barfmat)
write(mfile, "slackvals", slackmat)
close(mfile)






20 changes: 10 additions & 10 deletions codejl/input.jl
100644 → 100755
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ function readcase(casefilename, costsfilename, refbus, loadscale, thermalLimitSc
buses = Bus[]
for i in 1:size(busmat,1)
@assert busmat[i,1] == i # indexed in order
bustype = int(busmat[i,2])
bustype = round(Int,(busmat[i,2]))
Pd = busmat[i,3]
Qd = busmat[i,4]
Gs = busmat[i,5]
Expand All @@ -140,7 +140,7 @@ function readcase(casefilename, costsfilename, refbus, loadscale, thermalLimitSc
generatorlist = Int[]
genmat = case["gen"]
for i in 1:size(genmat,1)
busidx = int(genmat[i,1])
busidx = round(Int,genmat[i,1])
Pg = genmat[i,2]
Qg = genmat[i,3]
Pgmax = loadscale*genmat[i,9]
Expand All @@ -161,8 +161,8 @@ function readcase(casefilename, costsfilename, refbus, loadscale, thermalLimitSc
branchmat = case["branch"]
lines = Line[]
for i in 1:size(branchmat,1)
fbus = int(branchmat[i,1])
tbus = int(branchmat[i,2])
fbus = round(Int,branchmat[i,1])
tbus = round(Int,branchmat[i,2])
x = branchmat[i,4]
y = 1/x
u = branchmat[i,6]
Expand All @@ -181,7 +181,7 @@ function readcase(casefilename, costsfilename, refbus, loadscale, thermalLimitSc
if costsfilename != "none"
costsmat = readcsv(costsfilename, Float64)
for i in 1:size(costsmat,1)
busid = int(costsmat[i,1])
busid = round(Int,costsmat[i,1])
@assert length(buses[busid].genids) == 1
buses[busid].pi2 = costsmat[i,2] # quadratic coefficient
buses[busid].pi1 = costsmat[i,3] # linear coefficient
Expand All @@ -201,7 +201,7 @@ function readconfig(configfilename)
println("\nreading config $configfilename")
refbus = 0
uniformalphas = false

lines = readlines(open(configfilename,"r"))

numlines = length(lines)
Expand All @@ -221,8 +221,8 @@ function readconfig(configfilename)
extras = Dict()

for l in lines
beginswith(l,'#') && continue
startswith(l,'#') && continue

thisline = split(l)
length(thisline) > 0 || continue
if thisline[1] == "END"
Expand All @@ -234,7 +234,7 @@ function readconfig(configfilename)
elseif thisline[1] == "costs"
costsfilename = thisline[2]
elseif thisline[1] == "refbus"
refbus = int(thisline[2])
refbus = parse(Int,(thisline[2]))
elseif thisline[1] == "line_probability_threshold"
line_probability_threshold = float(thisline[2])
println(">>>> line_probability_threshold = $line_probability_threshold")
Expand Down Expand Up @@ -286,7 +286,7 @@ function readwind(windfilename, buses, loadscale)
push!(farms, Farm(int(id), loadscale*float(mean), loadscale*float(std)))
println("farm $i : $id $mean $std")
end

for i in 1:numfarms
setfarm(buses[farms[i].node],i)
end
Expand Down