Modeling with JuMP

This page shows a minimal example of PowerSystems.jl used to develop and Economic Dispatch model. The code shows the stages to develop modeling code

  1. Make the data set from power flow and time series data,
  2. Serialize the data,
  3. 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

The first step is to load the test data used throughout the rest of these tutorials (or set DATA_DIR as appropriate if it already exists).

julia> system_data =  build_system(PSISystems, "c_sys5_pjm")┌ Info: Building new system c_sys5_pjm from raw data
└   sys_descriptor.raw_data = "/home/runner/.julia/artifacts/2b7013812a985ce5884696ce276ca17b2719a57a/PowerSystemsTestData-2.0/psy_data/data_5bus_pu.jl"
[ Info: Serialized time series data to /home/runner/.julia/packages/PowerSystemCaseBuilder/8GyJc/data/serialized_system/ForecastOnly/c_sys5_pjm_time_series_storage.h5.
[ Info: Serialized System to /home/runner/.julia/packages/PowerSystemCaseBuilder/8GyJc/data/serialized_system/ForecastOnly/c_sys5_pjm.json
[ Info: Serialized System metadata to /home/runner/.julia/packages/PowerSystemCaseBuilder/8GyJc/data/serialized_system/ForecastOnly/c_sys5_pjm_metadata.json
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 │ Has Static Time Series │ Has Forecasts │
├───────────────────┼───────┼────────────────────────┼───────────────┤
│ ACBus             │ 5     │ false                  │ false         │
│ Arc               │ 6     │ false                  │ false         │
│ Line              │ 6     │ false                  │ false         │
│ PowerLoad         │ 3     │ true                   │ false         │
│ RenewableDispatch │ 2     │ true                   │ false         │
│ ThermalStandard   │ 5     │ false                  │ false         │
└───────────────────┴───────┴────────────────────────┴───────────────┘

Time Series Summary
┌──────────────────────────────────┬────────────┐
│ Property                         │ Value      │
├──────────────────────────────────┼────────────┤
│ Components with time series data │ 5          │
│ Total StaticTimeSeries           │ 5          │
│ Total Forecasts                  │ 0          │
│ Resolution                       │ 60 minutes │
└──────────────────────────────────┴────────────┘
julia> function ed_model(system::System, optimizer) 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) == net_load[t]) end @objective(ed_m, Min, sum( pg[get_name(g), t]^2*get_cost(get_variable(get_operation_cost(g)))[1] + pg[get_name(g), t]*get_cost(get_variable(get_operation_cost(g)))[2] for g in get_components(ThermalGen, system), t in time_periods ) ) optimize!(ed_m) return ed_m ended_model (generic function with 1 method)
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.14, running with linear solver MUMPS 5.6.2. 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.011 EXIT: Optimal Solution Found. A JuMP Model Minimization problem with: Variables: 120 Objective function type: JuMP.QuadExpr `JuMP.AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 24 constraints `JuMP.AffExpr`-in-`MathOptInterface.GreaterThan{Float64}`: 120 constraints `JuMP.AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 120 constraints `JuMP.VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 120 constraints Model mode: AUTOMATIC CachingOptimizer state: ATTACHED_OPTIMIZER Solver name: Ipopt Names registered in the model: pg