using Joulia
using DataFrames, CSV
using Gurobi
# data load for 2015 sample data
# see http://doi.org/10.5281/zenodo.1044463
pp_df = CSV.read("data_2015/power_plants.csv")
prices_df = CSV.read("data_2015/prices.csv")
# [...]
# generation of Joulia data types
pp = PowerPlants(pp_df, avail=avail_con_df, prices=prices_df)
storages = Storages(storages_df)
lines = Lines(lines_df)
nodes = Nodes(nodes_df, load_df, exchange_df)
res = RenewableEnergySource(res_df, avail)
# generation of the Joulia model
elmod = JouliaModel(pp, res, storages, nodes, lines)
# sclicing the data in weeks with the first full week starting at hour 49
slices = week_slices(49)
# running the Joulia model for week 30 using the Gurobi solver
results = run_model(elmod, slices[30], solver=GurobiSolver())
# src/model.jl
function create_model(pp::PowerPlants,
res::RenewableEnergySource,
nodes::Nodes)
N = nodes.id
P = pp.id
RES = res.id
mc = pp.mc
load = nodes.load
input_gen = pp.capacity
input_res = res.infeed
function build_model(T::Array{Int, 1})
m = Model()
@variables m begin
G[P,T] >= 0
G_RES[RES,T] >= 0
end
println("Set model objective")
@objective(m, Min, sum(mc[p][t] * G[p,t] for p in P, t in T))
println("Building constraints:")
prog = Progress(length(T), 0.2, "MarketClearing... ", 50)
@constraintref MarketClearing[T]
for t=T
MarketClearing[t] = @constraint(m,
sum(G[p,t] for p in P)
+ sum(G_RES[r,t] for r in RES)
==
sum(load[n][t] for n in N))
next!(prog)
end
JuMP.registercon(m, :MarketClearing, MarketClearing)
@constraintref GenerationRestriction[P,T]
prog = Progress(length(P)*length(T), 0.2, "Max generation... ", 50)
for p=P, t=T
GenerationRestriction[p,t] = @constraint(m,
G[p,t] <= input_gen[p][t] );
next!(prog)
end
@constraintref ResRestriction[RES,T]
prog = Progress(length(RES)*length(T), 0.2, "Renewable infeed... ", 50)
for r=RES, t=T
ResRestriction[r,t] = @constraint(m,
G_RES[r,t] <= input_res[r][t] );
next!(prog)
end
return m
end # end of build model
return build_model
end # end function `create_model`
println("Building constraints:")
prog = Progress(length(T), 0.2, "MarketClearing... ", 50)
@constraintref MarketClearing[T]
for t=T
MarketClearing[t] = @constraint(m,
sum(G[p,t] for p in P)
+ sum(G_RES[r,t] for r in RES)
==
sum(load[n][t] for n in N))
next!(prog)
end
JuMP.registercon(m, :MarketClearing, MarketClearing)
@constraintref GenerationRestriction[P,T]
prog = Progress(length(P)*length(T), 0.2, "Max generation... ", 50)
for p=P, t=T
GenerationRestriction[p,t] = @constraint(m,
G[p,t] <= input_gen[p][t] );
next!(prog)
end
run_opf("case3.m", DCPPowerModel, JuMP.with_optimizer(Gurobi.Optimizer))
is shorthand for
network_data = PowerModels.parse_file("case3.m")
pm = build_generic_model(network_data, ACPPowerModel, PowerModels.post_opf)
solve_generic_model(pm, JuMP.with_optimizer(Gurobi.Optimizer))
network_data holds the system parameters as stacked dictionaries (bus, gen, branch, ... as keys)ACPPowerModel = GenericPowerModel{ACPPowerForm} selects the formulation of the power flow to usepost_opf function defines the common recipe for building the model# src/prob/opf.jl
function post_opf(pm::GenericPowerModel)
variable_voltage(pm)
variable_generation(pm)
variable_branch_flow(pm)
variable_dcline_flow(pm)
objective_min_fuel_cost(pm)
constraint_voltage(pm)
for i in ids(pm, :ref_buses)
constraint_theta_ref(pm, i)
end
for i in ids(pm, :bus)
constraint_kcl_shunt(pm, i)
end
for i in ids(pm, :branch)
constraint_ohms_yt_from(pm, i)
constraint_ohms_yt_to(pm, i)
constraint_voltage_angle_difference(pm, i)
constraint_thermal_limit_from(pm, i)
constraint_thermal_limit_to(pm, i)
end
for i in ids(pm, :dcline)
constraint_dcline(pm, i)
end
end
# src/prob/opf.jl
function post_opf(pm::GenericPowerModel)
# [...]
variable_generation(pm)
# [...]
for i in ids(pm, :bus)
constraint_kcl_shunt(pm, i)
end
# [...]
end
# src/core/variable.jl
function variable_generation(pm::GenericPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
var(pm, nw, cnd)[:pg] = JuMP.@variable(pm.model,
[i in ids(pm, nw, :gen)],
lower_bound = ref(pm, nw, :gen, i, "pmin", cnd),
upper_bound = ref(pm, nw, :gen, i, "pmax", cnd)
)
# and the same for the reactive part
end
# src/form/dcp.jl
function constraint_kcl_shunt(pm::AbstractPowerModel{T}, i::Int, nw::Int, cnd::Int) where T <: AbstractDCPForm
bus_arcs = ref(pm, nw, :bus_arcs, i)
bus_gens = ref(pm, nw, :bus_gens, i)
bus_loads = ref(pm, nw, :bus_loads, i)
bus_pd = Dict(k => ref(pm, nw, :load, k, "pd", cnd) for k in bus_loads)
pg = var(pm, nw, cnd, :pg)
p = var(pm, nw, cnd, :p)
con(pm, nw, cnd, :kcl_p)[i] = JuMP.@constraint(pm.model,
sum(p[a] for a in bus_arcs) ==
sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)))
# omit reactive constraint
end
# src/form/acp.jl
function constraint_kcl_shunt(pm::GenericPowerModel{T}, i::Int, nw::Int, cnd::Int) where T <: AbstractACPForm
[...]
con(pm, n, c, :kcl_p)[i] = JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs))*vm^2)
con(pm, n, c, :kcl_q)[i] = JuMP.@constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs))*vm^2)
end
ps_dict = PowerSystems.parsestandardfiles("case5_re.m");
# ps_dict = PowerSystems.read_csv_data("<folder with gen.csv branch.csv bus.csv, .. load.csv>")
sys = PowerSystems.System(ps_dict)
System:
buses: Bus[Bus(name="bus1"), Bus(name="bus2"), Bus(name="bus3"), Bus(name="bus4"), Bus(name="bus5")]
generators:
GenClasses(T:5,R:2,H:0):
thermal: ThermalDispatch[ThermalDispatch(name="Solitude"), ThermalDispatch(name="Park City"), ThermalDispatch(name="Alta"), ThermalDispatch(name="Brighton"), ThermalDispatch(name="Sundance")]
renewable: RenewableCurtailment[RenewableCurtailment(name="SolarBusC"), RenewableCurtailment(name="WindBusA")]
hydro: nothing
(end generators)
loads: ElectricLoad[PowerLoad(name="bus2"), PowerLoad(name="bus3"), PowerLoad(name="bus4")]
branches: Branch[Line(name="1"), Line(name="2"), Line(name="3"), Line(name="4"), PhaseShiftingTransformer(name="5"), PhaseShiftingTransformer(name="6"), Line(name="7")]
storage: Storage[]
basepower: 100.0
time_periods: 25
sys.branches[1]
Line: name: 1 available: true connectionpoints: (from = Bus(name="bus1"), to = Bus(name="bus2")) r: 0.00281 x: 0.0281 b: (from = 0.00356, to = 0.00356) rate: 2.0 anglelimits: (min = -0.5235987755982988, max = 0.5235987755982988)
sys.generators.thermal[1]
ThermalDispatch: name: Solitude available: true bus: Bus(name="bus3") tech: TechThermal econ: EconThermal
sys.generators.thermal[1].tech
TechThermal: activepower: 3.24498 activepowerlimits: (min = 0.0, max = 5.2) reactivepower: 3.9 reactivepowerlimits: (min = -3.9, max = 3.9) ramplimits: (up = 0.0, down = 0.0) timelimits: nothing
sys.generators.thermal[1].econ
EconThermal:
capacity: 5.2
variablecost: getfield(PowerSystems, Symbol("##441#451")){Dict{String,Any}}(Dict{String,Any}("apf"=>0.0,"qc1max"=>0.0,"pg"=>3.24498,"model"=>2,"shutdown"=>0.0,"startup"=>0.0,"qc2max"=>0.0,"ramp_agc"=>0.0,"name"=>"Solitude","qg"=>3.9…))
fixedcost: 0.0
startupcost: 0.0
shutdncost: 0.0
annualcapacityfactor: nothing
head(sys.generators.renewable[1].scalingfactor)
6×1 TimeArray{Float64,1,DateTime,Array{Float64,1}} 2019-05-23T00:00:00 to 2019-05-23T05:00:00
│ │ A │
├─────────────────────┼───────┤
│ 2019-05-23T00:00:00 │ 1.0 │
│ 2019-05-23T01:00:00 │ 1.0 │
│ 2019-05-23T02:00:00 │ 1.0 │
│ 2019-05-23T03:00:00 │ 1.0 │
│ 2019-05-23T04:00:00 │ 1.0 │
│ 2019-05-23T05:00:00 │ 1.0 │
ps_dict = PowerSystems.parsestandardfiles("case5_re.m");
sys = PowerSystems.System(ps_dict);
econdispatch = PowerSimulations.EconomicDispatch(sys, PowerModels.DCPlosslessForm, optimizer=with_optimizer(Gurobi.Optimizer))
solve_op_model!(econdispatch)
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Warning, invalid warm-start basis discarded
Warning, invalid warm-start basis discarded
Optimize a model with 850 rows, 600 columns and 2025 nonzeros
Coefficient statistics:
Matrix range [1e+00, 2e+02]
Objective range [1e+03, 4e+03]
Bounds range [4e-01, 1e+01]
RHS range [5e-01, 4e+00]
Presolve removed 725 rows and 325 columns
Presolve time: 0.00s
Presolved: 125 rows, 275 columns, 550 nonzeros
Iteration Objective Primal Inf. Dual Inf. Time
0 0.0000000e+00 1.427467e+02 0.000000e+00 0s
191 4.3503066e+05 0.000000e+00 0.000000e+00 0s
Solved in 191 iterations and 0.00 seconds
Optimal objective 4.350306626e+05
PowerSimulations.OpertationModelResults(Dict(:Pre=>25×2 DataFrames.DataFrame │ Row │ SolarBusC │ WindBusA │ │ │ Float64 │ Float64 │ ├─────┼───────────┼──────────┤ │ 1 │ 0.6 │ 0.6 │ │ 2 │ 0.6 │ 0.6 │ │ 3 │ 0.6 │ 0.6 │ │ 4 │ 0.6 │ 0.6 │ │ 5 │ 0.6 │ 0.6 │ │ 6 │ 0.6 │ 0.6 │ │ 7 │ 0.6 │ 0.6 │ │ 8 │ 0.6 │ 0.6 │ │ 9 │ 0.6 │ 0.6 │ │ 10 │ 0.6 │ 0.6 │ ⋮ │ 15 │ 0.6 │ 0.6 │ │ 16 │ 0.6 │ 0.6 │ │ 17 │ 0.6 │ 0.6 │ │ 18 │ 0.6 │ 0.6 │ │ 19 │ 0.6 │ 0.6 │ │ 20 │ 0.6 │ 0.6 │ │ 21 │ 0.6 │ 0.6 │ │ 22 │ 0.6 │ 0.6 │ │ 23 │ 0.6 │ 0.6 │ │ 24 │ 0.6 │ 0.6 │ │ 25 │ 0.6 │ 0.6 │,:Pth=>25×5 DataFrames.DataFrame │ Row │ Solitude │ Park City │ Alta │ Brighton │ Sundance │ │ │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ ├─────┼──────────┼───────────┼─────────┼──────────┼──────────┤ │ 1 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 2 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 3 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 4 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 5 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 6 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 7 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 8 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 9 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 10 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ ⋮ │ 15 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 16 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 17 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 18 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 19 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 20 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 21 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 22 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 23 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 24 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │ │ 25 │ 3.8 │ 1.7 │ 0.4 │ 2.9 │ 0.0 │), Dict(:ED=>4.35031e5), Dict{Symbol,Any}(:dual_status=>FEASIBLE_POINT,:primal_status=>FEASIBLE_POINT,:termination_status=>OPTIMAL))
devices = Dict{Symbol, PSI.DeviceModel}(:ThermalGenerators => PSI.DeviceModel(PSY.ThermalGen, PSI.ThermalDispatch),
:RenewableGenerators => PSI.DeviceModel(PSY.RenewableGen, PSI.RenewableFullDispatch),
:Loads => PSI.DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))
branches = Dict{Symbol, PSI.DeviceModel}(:Lines => PSI.DeviceModel(PSY.Branch, PSI.SeriesLine))
# [...]
econdispatch = PowerOperationModel(EconomicDispatch,
PowerModels.DCPlosslessForm,
devices,
branches,
services,
system,
optimizer = with_optimizer(Gurobi.Optimizer))
fieldnames(typeof(econdispatch.canonical_model))
(:JuMPmodel, :variables, :constraints, :cost_function, :expressions, :parameters, :initial_conditions, :pm_model)
# src/devices/device_constructors/thermalgeneration_constructor.jl
function construct_device!(ps_m::CanonicalModel,
device::Type{T},
device_formulation::Type{D},
system_formulation::Type{S},
sys::PSY.System,
time_range::UnitRange{Int64};
kwargs...) where {T<: PSY.ThermalGen,
D <: AbstractThermalDispatchForm,
S <: PM.AbstractActivePowerFormulation}
#Variables
activepower_variables(ps_m, sys.generators.thermal, time_range);
#Constraints
activepower_constraints(ps_m, sys.generators.thermal, device_formulation, system_formulation, time_range)
#Cost Function
cost_function(ps_m, sys.generators.thermal, device_formulation, system_formulation)
return
end
# in src/devices/device_models/thermal_generation.jl
function activepower_constraints(ps_m::CanonicalModel,
devices::Array{T,1},
device_formulation::Type{D},
system_formulation::Type{S},
time_range::UnitRange{Int64}) where {T <: PSY.ThermalGen,
D <: AbstractThermalDispatchForm,
S <: PM.AbstractPowerFormulation}
range_data = [(g.name, g.tech.activepowerlimits) for g in devices]
device_range(ps_m, range_data, time_range, :thermal_active_range, :Pth)
end
# in src/devices/device_models/common/range_constraint.jl -- simplified
function device_range(ps_m::CanonicalModel,
range_data::Vector{NamedMinMax},
time_range::UnitRange{Int64},
cons_name::Symbol,
var_name::Symbol)
ps_m.constraints[cons_name] = JuMP.@constraint(ps_m.JuMPmodel, [r=range_data,t=time_range], r.min <= ps_m.variables[var_name][r, t] <= r.max)
end
model = EnergyModel("elec_s_45_lv1.25_3H.nc", solver=GurobiSolver())
build(model);
solve(model)
Everything is fetched on demand from NetCDF
get(model, :onwind, :p_nom) # -> AxisArray of data in netcdf file variable onwind::p_nom
Component
├──OnePort
│ ├──Generator
│ ├──StorageUnit
│ └──Store
├──PassiveBranch
│ ├──Line
│ └──Transformer
└──ActiveBranch
└──Link# src/components/oneport.jl
cost(c::OnePort) = sum(c[:marginal_cost] .* c[:p]) + sum(c[:capital_cost] .* (c[:p_nom] - getparam(c, :p_nom)))
function nodalbalance(c::OnePort)
p = c[:p]
c[:bus] => (o,t)->p[o,t]
end
## Generator
function build(c::Generator)
T = axis(c, :snapshots)
G = axis(c)
if isvar(c, :p_nom)
@emvariable c c[:p_nom_min][g] <= p_nom[g=G] <= c[:p_nom_max][g]
end
@emvariable c c[:p_min_pu][g,t] * c[:p_nom][g] <= p[g=G,t=T] <= c[:p_max_pu][g,t] * c[:p_nom][g]
end
addelement(Generator, :generators, (:G, :T=>:snapshots), joinpath(@__DIR__, "generators.csv"))
# src/components/generators.csvcsv
attribute,quantitytype,dimensions,default,dtype
p_nom,VarParam,G,0,float
p_max_pu,Param,"G,T",1,float
p_min_pu,Param,"G,T",0,float
p_nom_min,Param,G,0,float
p_nom_max,Param,G,Inf,float
efficiency,Param,G,1,float
carrier,Param,,,string
marginal_cost,Param,"G,T",0.,float
capital_cost,Param,G,0.,float
bus,Param,G,,string
p,Variable,"G,T",0,float| Joulia | PowerModels | PowerSystems/PowerSimulation | EnergyModels | |
|---|---|---|---|---|
| Storage of | ||||
| - Parameters | Flat types: Nodes, Lines, PowerPlants, RenewableEnergySource, Storages, which contain dictionaries for each parameter indexed by names |
Matpower files are parsed into a nested structure of dictionaries. ref() - Getter helps with access |
Nested detailed type hierarchy to hold data for individual components in PowerSystems.jl, possibility for extension, importers from matpower and raw | Parameters are fetched on demand from a NetCDF backend as flexibly dimensioned DenseAxisArrays |
| - Variables | Inside JuMP Model: model[:MarketClearing] |
Variables and Constraints are equally sorted into nested dictonaries using the helper functions var() and con() |
Stored in dictionary on canonical_model.{variables,constraints} mapping a chosen symbols to DenseAxisArrays |
Variables and constraint are stored on the JuMP model with a dedicated prefix involving the class name |
| Primary focus/features | Economic dispatch with thermal and renewable generators, storage units including linearized powerflows (based on PTDF) | Framework for reference implementations and easy comparison of several power flow formulations, with extensions for several timesteps, storage und network expansion | Production cost modelling with flexible device model representations, Priority focus has been on operational models, Capacity expansion is coming up, integration with PowerDynamics is planned | Co-Optimization of transmission, storage and generation capacities under linearized optimal power flow constraints for large networks |
install julia by getting it from http://julialang.org/ and extracting it to a directory of your choice.
For windows during the workshop there was a comment that you have to adjust your PATH environment variable to include Julia's bin directory for IJulia to work.
On linux, I tend to unpack Julia to ~/.julia/julia-1.1 for instance, and then symlink the julia binary to ~/.local/bin (ie. ln -s ~/.julia/julia-1.1/bin/julia ~/.local/bin/julia)
After installation you can launch Julia (by clicking or calling julia) and add its kernel with
]add IJulia
(the closing bracket switches into Julia's package manager mode)
There is an example to run elmod in Joulia at https://github.com/JuliaEnergy/JouliaExamples.
For now, you will have to install Joulia using:
]add JuMP@v0.18
]add https://github.com/JuliaEnergy/Joulia.jlNREL compiled notebooks for demonstrating their libraries in the repository: https://github.com/NREL-SIIP/Examples
Once you got/cloned the repository and changed into it, you can get into some sort an environment where everything is installed with
]activate env; instantiateEnergyModels needs an overhaul to work with the current Julia versions. So to use it at the moment you need to install Julia 0.6. An update is coming soon.
using Pkg
Pkg.add("https://gitlab.com/coroa/EnergyModels.git")
]add PowerModelsThere are test systems for PowerModels and PowerSystems in the repository
https://github.com/NREL/PowerSystemsTestData
PowerModels reads primarily the matpower format.