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:
- 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.
Start by loading required packages:
julia> using PowerSystemsjulia> using JuMPjulia> using Ipoptjulia> using PowerSystemCaseBuilderjulia> using Dates
For this example, we'll load an existing data set using PowerSystemCaseBuilder.jl, which is a helper library that makes it easier to reproduce examples. Normally you would pass your local files to create the system data instead of calling the function build_system. We also use transform_single_time_series! to format time-series data as forecasts for this problem:
julia> system_data = build_system(PSISystems, "c_sys5_pjm")[ Info: Loaded time series from storage file existing=/home/runner/.julia/packages/PowerSystemCaseBuilder/QOo1f/data/serialized_system/614e094ea985a55912fc1321256a49b755f9fc451c0f264f24d6d04af20e84d7/c_sys5_pjm_time_series_storage.h5 new=/tmp/jl_xL7kxT compression=CompressionSettings(false, CompressionTypes.DEFLATE = 1, 3, true) 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 │ name │ time_series_type │ i ⋯ │ String │ String │ String │ String │ S ⋯ ├───────────────────┼────────────────┼──────────────────┼──────────────────┼──── │ PowerLoad │ Component │ max_active_power │ SingleTimeSeries │ 2 ⋯ │ RenewableDispatch │ Component │ max_active_power │ SingleTimeSeries │ 2 ⋯ └───────────────────┴────────────────┴──────────────────┴──────────────────┴──── 4 columns omittedjulia> transform_single_time_series!(system_data, Hour(24), Hour(24))
Next, we define the custom 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> 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 ended_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.19, running with linear solver MUMPS 5.8.1. 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.008 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