Modeling with JuMP
This guide is for users who are interested in writing custom optimization problems directly in JuMP, using data formatted with PowerSystems.jl
. Check out PowerSimulations.jl
for developing reusable templates for optimization problems within the Sienna platform.
This page shows a minimal example to develop a Economic Dispatch model. The code shows the stages to develop modeling code:
This page shows a minimal example to develop a Economic Dispatch model. The code shows the stages to develop modeling code:
- Make the data set from power flow and time series data,
- Serialize the system data,
- Pass the data and algorithm to the model.
One of the main uses of PowerSystems.jl
is not having re-run the data generation for every model execution. The model code shows an example of populating the constraints and cost functions using accessor functions inside the model function. The example concludes by reading the data created earlier and passing the algorithm with the data.
julia> using PowerSystems
julia> const PSY = PowerSystems
PowerSystems
julia> using JuMP
julia> using Ipopt
julia> using PowerSystemCaseBuilder
julia> using Dates
PowerSystemCaseBuilder.jl
is a helper library that makes it easier to reproduce examples in the documentation and tutorials. Normally you would pass your local files to create the system data instead of calling the function build_system
. For more details visit PowerSystemCaseBuilder Documentation.
Next, we define the optimization problem using JuMP
's syntax. The constraints include each generator's minimum and maximum active power output as well as the system power balance equation, minimizing the operating cost for each step in the 24-hour horizon:
julia> system_data = build_system(PSISystems, "c_sys5_pjm")
[ Info: Loaded time series from storage file existing=/home/runner/.julia/packages/PowerSystemCaseBuilder/dJGb8/data/serialized_system/614e094ea985a55912fc1321256a49b755f9fc451c0f264f24d6d04af20e84d7/c_sys5_pjm_time_series_storage.h5 new=/tmp/jl_ITTqoY compression=CompressionSettings(false, CompressionTypes.DEFLATE = 1, 3, true) [ Info: rating 1114.8 MW for 4 is outside the expected range (min = 327.0, max = 797.0) MW for Line at a 230.0 kV Voltage level. [ Info: rating 200.0 MW for 1 is outside the expected range (min = 327.0, max = 797.0) MW for Line at a 230.0 kV Voltage level. ┌ Warning: rating 4053.0 MW for 5 is 2x larger than the max expected rating 797.0 MW for Line at a 230.0 kV Voltage level. └ @ PowerSystems ~/work/PowerSystems.jl/PowerSystems.jl/src/utils/IO/branchdata_checks.jl:102 [ Info: rating 200.0 MW for 2 is outside the expected range (min = 327.0, max = 797.0) MW for Line at a 230.0 kV Voltage level. [ Info: rating 200.0 MW for 6 is outside the expected range (min = 327.0, max = 797.0) MW for Line at a 230.0 kV Voltage level. ┌ Warning: rating 1881.2 MW for 3 is 2x larger than the max expected rating 797.0 MW for Line at a 230.0 kV Voltage level. └ @ PowerSystems ~/work/PowerSystems.jl/PowerSystems.jl/src/utils/IO/branchdata_checks.jl:102 System ┌───────────────────┬─────────────┐ │ Property │ Value │ ├───────────────────┼─────────────┤ │ Name │ │ │ Description │ │ │ System Units Base │ SYSTEM_BASE │ │ Base Power │ 100.0 │ │ Base Frequency │ 60.0 │ │ Num Components │ 27 │ └───────────────────┴─────────────┘ Static Components ┌───────────────────┬───────┐ │ Type │ Count │ ├───────────────────┼───────┤ │ ACBus │ 5 │ │ Arc │ 6 │ │ Line │ 6 │ │ PowerLoad │ 3 │ │ RenewableDispatch │ 2 │ │ ThermalStandard │ 5 │ └───────────────────┴───────┘ StaticTimeSeries Summary ┌───────────────────┬────────────────┬──────────────────┬─────────────────────── │ owner_type │ owner_category │ time_series_type │ initial_timestamp ⋯ │ String │ String │ String │ String ⋯ ├───────────────────┼────────────────┼──────────────────┼─────────────────────── │ PowerLoad │ Component │ SingleTimeSeries │ 2024-01-01T00:00:00 ⋯ │ RenewableDispatch │ Component │ SingleTimeSeries │ 2024-01-01T00:00:00 ⋯ └───────────────────┴────────────────┴──────────────────┴─────────────────────── 3 columns omitted
julia> transform_single_time_series!(system_data, Hour(24), Hour(24))
julia> function ed_model(system::System, optimizer, load_scaling_factor::Float64 = 1.0) ed_m = Model(optimizer) time_periods = 1:24 thermal_gens_names = get_name.(get_components(ThermalStandard, system)) @variable(ed_m, pg[g in thermal_gens_names, t in time_periods] >= 0) for g in get_components(ThermalStandard, system), t in time_periods name = get_name(g) @constraint(ed_m, pg[name, t] >= get_active_power_limits(g).min) @constraint(ed_m, pg[name, t] <= get_active_power_limits(g).max) end net_load = zeros(length(time_periods)) for g in get_components(RenewableGen, system) net_load -= get_time_series_values(SingleTimeSeries, g, "max_active_power")[time_periods] end for g in get_components(StaticLoad, system) net_load += get_time_series_values(SingleTimeSeries, g, "max_active_power")[time_periods] end for t in time_periods @constraint( ed_m, sum(pg[g, t] for g in thermal_gens_names) == load_scaling_factor * net_load[t] ) end @objective( ed_m, Min, sum( pg[get_name(g), t] * get_proportional_term(get_function_data(get_variable(get_operation_cost(g)))) for g in get_components(ThermalGen, system), t in time_periods ) ) optimize!(ed_m) return ed_m end
ed_model (generic function with 2 methods)
Finally, the PowerSystems.jl
data is combined with this economic dispatch model and solved with the open-source Ipopt
solver:
julia> results = ed_model(system_data, Ipopt.Optimizer)
****************************************************************************** This program contains Ipopt, a library for large-scale nonlinear optimization. Ipopt is released as open source code under the Eclipse Public License (EPL). For more information visit https://github.com/coin-or/Ipopt ****************************************************************************** This is Ipopt version 3.14.17, running with linear solver MUMPS 5.7.3. Number of nonzeros in equality constraint Jacobian...: 120 Number of nonzeros in inequality constraint Jacobian.: 240 Number of nonzeros in Lagrangian Hessian.............: 0 Total number of variables............................: 120 variables with only lower bounds: 120 variables with lower and upper bounds: 0 variables with only upper bounds: 0 Total number of equality constraints.................: 24 Total number of inequality constraints...............: 240 inequality constraints with only lower bounds: 120 inequality constraints with lower and upper bounds: 0 inequality constraints with only upper bounds: 120 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 0 2.6159974e+01 2.52e+00 6.07e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 6.4569853e+02 7.10e-01 3.42e+01 -1.0 5.65e-01 - 2.08e-02 7.18e-01h 1 2 5.2023163e+02 6.87e-01 3.31e+01 -1.0 2.09e+01 - 2.97e-02 3.25e-02f 1 3 4.4033326e+02 6.50e-01 3.14e+01 -1.0 6.64e+00 - 1.08e-01 5.30e-02f 1 4 4.0560222e+02 5.78e-01 2.79e+01 -1.0 2.88e+00 - 1.73e-01 1.12e-01f 1 5 4.0100081e+02 5.39e-01 2.60e+01 -1.0 5.67e+00 - 9.29e-02 6.75e-02f 1 6 3.9908466e+02 4.48e-01 2.16e+01 -1.0 1.99e+00 - 3.84e-01 1.68e-01h 1 7 4.3287862e+02 1.43e-01 7.23e+00 -1.0 5.96e-01 - 1.00e+00 6.81e-01h 1 8 4.6122880e+02 4.44e-16 1.00e-06 -1.0 1.51e-01 - 1.00e+00 1.00e+00h 1 9 4.4374486e+02 4.44e-16 2.63e-02 -2.5 9.43e-02 - 9.17e-01 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 10 4.4264070e+02 4.44e-16 1.50e-09 -3.8 9.08e-03 - 1.00e+00 1.00e+00f 1 11 4.4261103e+02 4.44e-16 1.85e-11 -5.7 1.71e-04 - 1.00e+00 1.00e+00f 1 12 4.4261068e+02 4.44e-16 2.84e-14 -8.6 1.98e-06 - 1.00e+00 1.00e+00f 1 Number of Iterations....: 12 (scaled) (unscaled) Objective...............: 4.4261067891368759e+02 4.4261067891368759e+02 Dual infeasibility......: 2.8421709430404007e-14 2.8421709430404007e-14 Constraint violation....: 4.4408920985006262e-16 4.4408920985006262e-16 Variable bound violation: 9.8329323571942150e-09 9.8329323571942150e-09 Complementarity.........: 2.5091354477996078e-09 2.5091354477996078e-09 Overall NLP error.......: 2.5091354477996078e-09 2.5091354477996078e-09 Number of objective function evaluations = 13 Number of objective gradient evaluations = 13 Number of equality constraint evaluations = 13 Number of inequality constraint evaluations = 13 Number of equality constraint Jacobian evaluations = 1 Number of inequality constraint Jacobian evaluations = 1 Number of Lagrangian Hessian evaluations = 1 Total seconds in IPOPT = 0.114 EXIT: Optimal Solution Found. A JuMP Model ├ solver: Ipopt ├ objective_sense: MIN_SENSE │ └ objective_function_type: JuMP.AffExpr ├ num_variables: 120 ├ num_constraints: 384 │ ├ JuMP.AffExpr in MOI.EqualTo{Float64}: 24 │ ├ JuMP.AffExpr in MOI.GreaterThan{Float64}: 120 │ ├ JuMP.AffExpr in MOI.LessThan{Float64}: 120 │ └ JuMP.VariableRef in MOI.GreaterThan{Float64}: 120 └ Names registered in the model └ :pg