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.

julia> using PowerSystems
julia> const PSY = PowerSystemsPowerSystems
julia> using JuMP
julia> using Ipopt
julia> using PowerSystemCaseBuilder
julia> using Dates

Note

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
       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.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