# Hierarchical/compositional construction of large optimization models

**URL:** <https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303>\
**Category:** Optimization (Mathematical)\
**Created:** [June 13, 2023, 9:30pm UTC](https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303 "2023-06-13T21:30:34Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![slwu89](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/slwu89/32/217323_2.png) [@slwu89](https://discourse.julialang.org/u/slwu89)\
**Post date:** [June 13, 2023, 9:30pm UTC](https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303/1 "2023-06-13T21:30:34Z")

</div>

Hi everyone. I have a broad question, rather than a question about any specific package or code block, so please let me know if its inappropriate and I’ll move/close it.

I’m interested if the Julia/JuMP community has published/learned any best practices for setting up large models for optimization with a naturally hierarchical structure. I’ve searched around for examples to learn best practices from and came up with a somewhat overwhelming list of packages which set up large JuMP models, all having to do with power infrastructure (the specific problems I am interested in are not related to power generation/distribution, otherwise I’d just use one of these):

- [leonardgoeke/AnyMOD.jl: Julia framework for energy system models with a focus on multi-period capacity expansion (github.com)](https://github.com/leonardgoeke/AnyMOD.jl)
- [JuliaEnergy/PowerDynamics.jl: Package for dynamical modeling of power grids (github.com)](https://github.com/JuliaEnergy/PowerDynamics.jl)
- [Breakthrough-Energy/REISE.jl: Renewable Energy Integration Simulation Engine (github.com)](https://github.com/Breakthrough-Energy/REISE.jl)
- [NREL-Sienna/PowerSystems.jl: Data structures in Julia to enable power systems analysis. Part of the Scalable Integrated Infrastructure Planning Initiative at the National Renewable Energy Lab. (github.com)](https://github.com/NREL-Sienna/PowerSystems.jl)
- [lanl-ansi/PowerModels.jl: A Julia/JuMP Package for Power Network Optimization (github.com)](https://github.com/lanl-ansi/PowerModels.jl)
- [andrewrosemberg/HydroPowerModels.jl: A Julia/JuMP Package for Hydrothermal economic dispatch Optimization (github.com)](https://github.com/andrewrosemberg/HydroPowerModels.jl)
- [martinbiel/HydroModels.jl: Julia modeling framework for hydropower operations. (github.com)](https://github.com/martinbiel/HydroModels.jl)
- [psrenergy/SDDPDSO.jl: SDDP-DSO stands for Stochastic Dual Dynamic Programming for Distribution System Operation (github.com)](https://github.com/psrenergy/SDDPDSO.jl/tree/master)
- [plasmo-dev/Plasmo.jl: A Platform for Scalable Modeling and Optimization (github.com)](https://github.com/plasmo-dev/Plasmo.jl)
- [GenXProject/GenX: GenX: a configurable power system capacity expansion model for studying low-carbon energy futures. More details at : https://genx.mit.edu (github.com)](https://github.com/GenXProject/GenX)
- [NREL/REopt.jl (github.com)](https://github.com/NREL/REopt.jl)

Are there any other packages I missed that the community would recommend as a good example of “best practices” to build large hierarchical problems in JuMP? And are any of the examples here exceptionally worth taking a look at the design to incorporate best practices? Are there any other references I should be checking out related to setting up large problems like this (I’ve read the “design patterns for larger models” page of the JuMP manual).

I realize I mixed a few packages that set up stochastic problems, relying on SDDP. In general I’m interested in both approaches, although mainly starting with the deterministic setting to begin with.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 13, 2023, 9:59pm UTC](https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303/2 "2023-06-13T21:59:38Z")

</div>

For other readers, here’s the link to the JuMP docs, which summarizes my suggested structure:

[Design patterns for larger models · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/getting_started/design_patterns_for_larger_models/)?

That list is fairly comprehensive. The NREL/SIIP PowerSystems.jl is probably the largest example with the most engineering time put into it. They went full modularity, and full performance (sometimes at the expense of some complicated code to pre-allocate JuMP containers).

The power-focused applications are also somewhat niche, because they have some unique features that aren’t in other applications (like a transmission graph, or multiple time structures, and very flexible generator combinations with different technologies all producing the same commodity).

I think the most important thing to start with is an understanding of the data. Ultimately, the JuMP optimization part will be a very small component of the overall application.

- What at the the inputs and outputs at each layer of the hierarchy?
- What format should the inputs and outputs be?
- How can you setup, test, and validate each layer of the hierarchy independently?

I’d also advise that unless necessary, you avoid fancy tricks. Do the simplest thing possible, even if it is slightly slower than a more clever thing.

As a cautionary tale of something that works, but which I don’t think we got right (and since I’m somewhat responsible), take a look at [GitHub - EPOC-NZ/JADE.jl](https://github.com/EPOC-NZ/JADE.jl).

- The inputs and outputs are just CSV files: [JADE.jl/test/Input/test1 at main · EPOC-NZ/JADE.jl · GitHub](https://github.com/EPOC-NZ/JADE.jl/tree/main/test/Input/test1)
- but they aren’t actually CSV files because various bits got added to them: [JADE.jl/test/Input/test1/thermal\_fuel\_costs.csv at main · EPOC-NZ/JADE.jl · GitHub](https://github.com/EPOC-NZ/JADE.jl/blob/main/test/Input/test1/thermal_fuel_costs.csv)
- We can’t test each piece in isolation, we need to run a full optimize/simulate workflow to check a single objective value: [JADE.jl/test/test\_cases.jl at f2e753eac0f6d7354c3b412da830c0d4c1bcec9c · EPOC-NZ/JADE.jl · GitHub](https://github.com/EPOC-NZ/JADE.jl/blob/f2e753eac0f6d7354c3b412da830c0d4c1bcec9c/test/test_cases.jl#L44-L56)
- There are lots of nested `if` statements as new features got requested: [JADE.jl/src/simulate.jl at f2e753eac0f6d7354c3b412da830c0d4c1bcec9c · EPOC-NZ/JADE.jl · GitHub](https://github.com/EPOC-NZ/JADE.jl/blob/f2e753eac0f6d7354c3b412da830c0d4c1bcec9c/src/simulate.jl#L150-L187)

Despite all the mess, the actual SDDP model is 600 lines, but surprisingly clean and easy to read and understand: [JADE.jl/src/model.jl at f2e753eac0f6d7354c3b412da830c0d4c1bcec9c · EPOC-NZ/JADE.jl · GitHub](https://github.com/EPOC-NZ/JADE.jl/blob/f2e753eac0f6d7354c3b412da830c0d4c1bcec9c/src/model.jl#L56-L620)

The backstory is the Julia code was a 2016 re-implementation of a C++ application from 2011, which was a re-implementation of a series of AMPL scripts from 2008. And then between 2016 and now, a bunch of things got added to it in various places. (And during that time, Julia and JuMP and SDDP.jl all changed quite a lot.). So at no point did we ever step back and design a nice representation of the data layer.

I guess my point is ignore the JuMP part to start with. If you have good data representation, the JuMP part will be clean and simple. If you have a nice JuMP model and messy data, the entire thing will be a mess.

---

<div class="post-metadata">

**Author:** ![slwu89](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/slwu89/32/217323_2.png) [@slwu89](https://discourse.julialang.org/u/slwu89)\
**Post date:** [June 13, 2023, 10:40pm UTC](https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303/3 "2023-06-13T22:40:07Z")

</div>

Thanks very much for the great advice!

With regards to the point of not being able to test each piece in isolation in JADE.jl, would you have preferred to set up the model/data input such that you could create SDDP/JuMP models for “parts” of the overall model and test them individually regarding constraints that link them to other bits of the overall model as input? For my problem, I’m really interested in being able to test in this way.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 13, 2023, 11:01pm UTC](https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303/4 "2023-06-13T23:01:43Z")

</div>

> would you have preferred to set up the model/data input such that you could create SDDP/JuMP models for “parts” of the overall model and test them individually regarding constraints that link them to other bits of the overall model as input?

💯 but I don’t know if I have a good answer for exactly what that looks like. It’s something that I’d like to address once the nonlinear rewrite is finished.

Potentially something like:

```julia
using JuMP

abstract type AbstractConstraint end

struct Data
    constraint_types::Vector{AbstractConstraint}
end

function initialize_model(data)
    return Model()
end

function add_variables(model, data)
    @variable(model, x)
end

struct ConstraintA <: AbstractConstraint end

function add_constraint(model, data, ::ConstraintA)
    @constraint(model, model[:x] <= 1)
end

function test_ConstraintA()
    data = Data([ConstraintA()])
    model = initialize_model(data)
    add_variables(model, data)
    add_constraint(model, data, ConstraintA())
    # https://jump.dev/JuMP.jl/dev/reference/solutions/#JuMP.primal_feasibility_report
    @test primal_feasibility_report(model, Dict(model[:x] => 0.0)) === nothing
    @test primal_feasibility_report(model, Dict(model[:x] => 1.0)) === nothing
    @test primal_feasibility_report(model, Dict(model[:x] => 1.01)) !== nothing
    @test primal_feasibility_report(model, Dict(model[:x] => 2.0)) !== nothing
    return
end

function main(data)
    model = initialize_model(data)
    add_variables(model, data)
    for constraint in data.constraint_types
        add_constraint(model, data, constraint)
    end
    return model
end

```

But ideally simpler, so that you build the entire model and then pass good/pass points for each constraint.

For smaller applications, this might be too much complication. Everything’s a trade-off. Getting the input data into the right format and validating it is more important than the JuMP details though.

---

<div class="post-metadata">

**Author:** ![slwu89](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/slwu89/32/217323_2.png) [@slwu89](https://discourse.julialang.org/u/slwu89)\
**Post date:** [September 4, 2024, 4:20am UTC](https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303/5 "2024-09-04T04:20:56Z")

</div>

@odow I hope you don’t mind me bumping this old conversation. I was curious if after some time there were any developments regarding your thoughts on the modular construction of things (especially since the nonlinear rewrite ended up finished).

After a year of experience in this (at that time) new field to me, I completely agree that the data modeling stage is often more important than the details of building the optimization model.

That being said, as you may have noticed from my response in another thread (linked below), I’m very interested in the potential of the C-Set data structure to bring the data model and the optimization model closer together. Hence my interest to understand in a bit more detail the thoughts you sketched out in your last reply here, I’m interested to see if this approach can also have benefits for this modular construction/testing of the optimization model.

> [@JuMP takes long time to load the problem/model](https://discourse.julialang.org/t/jump-takes-long-time-to-load-the-problem-model/118593/14):
>
> Very nice new tutorial. Being interested in model building for these things, I typed up a quick tutorial of using acsets (attributed C-sets, a data structure from applied category theory, but in practice something akin to an in-memory relational database with additional equations on paths possible) for generating this type of model. They are exported from the Catlab.jl package. The full script is here [https://jump.dev/JuMP.jl/dev/tutorials/getting\_started/sum\_if/ with acsets (github.com)](https://gist.github.com/slwu89/6adf107fb345a31101545520362789e7) The da…

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [September 4, 2024, 5:34am UTC](https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303/6 "2024-09-04T05:34:47Z")

</div>

I think the design patterns tutorial is still the best reference.

Re the nonlinear: there wasn’t anything specific to that that i had in mind. Just that i was busy with that development and didn’t have time for other things.

There are two open issues you might be interested in

> <https://github.com/jump-dev/JuMP.jl/issues/3664>
>
> Despite making it easy for users to formulate and solve a wide range of optimiza…tion problems, JuMP provides little support for the users who make mistakes, or tools for advanced users to debug problematic models. Moreover, in our experience the majority of (expensive) human programmer time is spent, not in formulating or solving a model, but in the debugging and testing stage of development.
> 
> Common questions that developers ask us during the debugging and testing stage of development include:
> 
> \* How do I find out why an optimization model is infeasible or unbounded?
> \* How do I know if the model has numerical issues?
> \* How to identify parts of the model that could be improved?
> 
> The goal of this project is to develop an automated analysis framework that can identify common problems to reduce the time spent by developers in the debugging and testing phase. There are a range of problems which are trivial to identify and fix. We suggest that the correct approach is not to fix the problems automatically (as the AMPL presolve does), but to educate the user on how to write better models.
> 
> \## Ideas
> 
> The deliverable for this issue is a Julia package that exposes an API for linting JuMP and MathOptInterface models. See below for a set of initial ideas.
> 
> \### Summary statistics
> 
> Print summary statistics of the model such as the number of variables and constraints, whether bridging is used, and a hash of the problem data.
> 
> Motivation: Simple statistics are often helpful to give an idea of the model. If the number of variables or constraints is very large, it may be difficult to solve, and if the number of variables or constraints does not match what the user expected, they may have made a mistake in the formulation.
> 
> Prior art: Most solvers provide this information. It is trivial to implement.
> 
> \### Coefficient analysis
> 
> For each constraint type, print statistics such as the number of constraints, the range of left- and right-hand side coefficients, and the density of the coefficients (number of terms / number of variables in the model), and whether the function is convex or non-convex.
> 
> Motivation: Poorly scaled problems have coefficients that vary by many orders of magnitude. Poor problem scaling is the most common cause of numerical problems as many common solution algorithms include a decomposition of the constraint matrix whose stability depends on the condition number of that matrix. Poorly scaled problems are usually an indication of user-error, for example, using the wrong units so a coefficient is off by a factor of 10^N. Once problematic coefficients are identified, they can usually be fixed by scaling the variables and input coefficients.
> 
> Prior art: Some solvers provide this information. Oscar has implemented it before in SDDP.jl: https://github.com/odow/SDDP.jl/blob/9e88c6ba2a7dacc0f842b9be3f8073acfd55073e/src/print.jl#L191-L384 
> 
> \### Degeneracy
> 
> Identify constraints in the problem that are linearly dependent. 
> 
> Motivation: Most solution algorithms require a linearly independent set of constraints. Because linearly dependent constraints are commonly provided by users, most solvers have detection routines for them. However, they are usually a sign that a model’s formulation can be improved.
> 
> Prior art: Alex Dowling has a paper and implementation https://github.com/adowling2/DegeneracyHunter.jl
> 
> \### Incidence analysis
> 
> Identify nonlinear models in which the Jacobian is not full rank. This problem is closely related to the degeneracy problem.
> 
> Motivation: Under- or over-constrained systems of nonlinear equations cause solvers trouble, and they are typically a sign of user-error since physical systems are typically not over- or under-constrained.
> 
> Prior art: Robby Parker has an implementation of this for Pyomo https://pyomo.readthedocs.io/en/stable/contributed\_packages/incidence/index.html, and an in-development version of the code for JuMP.
> 
> \### Bounds given as constraints
> 
> Identify constraints of the form \`@constraint(model, l \<= x \<= u)\`. These should instead be specified as \`@variable(model, l \<= x \<= u)\`. 
> 
> Motivation: The former adds a new linear constraint row to the problem. The latter sets a variable’s bounds. Solvers typically have specialised support for variable bounds. If the solver doesn’t have specialised support for variable bounds, JuMP will reformulate into the constraint version.
> 
> Prior art: Most solvers have a presolve that does this, but it’s a sign of user-error that could be improved.
> Variables that do not appear in constraints
> Variables that do not appear in constraints are generally a sign of user-error. They might have forgotten to include the variable, or they could make the model more performant by omitting the variable in the first place.
> 
> \### Starting point analysis
> 
> Starting points which violate domain restrictions like \`log(x)\` when \`x = 0\` to start are common. Therefore, we should analyse the feasibility of a given starting point, and potentially that of the first- and second-derivatives as well.
> 
> \### Domain analysis
> 
> Nonlinear programs often include terms like \`log(x)\` where \`x \>= 0\`. This is a common cause of domain violations because the solver may try \`x = 0\`. It’s usually better to add an explicit lower bound on \`x\` of \`x \>= 0.01\`.
> 
> \### Convexity analysis 
> 
> Proving convexity is difficult. But we could test feasible trial points for a proof of nonconvexity. A full disciplined convex programming package is out-of-scope for this project.
> 
> Prior art: CVXPY has an implementation: https://github.com/cvxpy/cvxpyanalyzer/blob/master/analyzer/convexity\_checker.py
> 
> \### Irreducible infeasible subsystem
> 
> If the problem is infeasible we could try and find a reduced subsystem that is still infeasible. This does not need to be the minimal subsystem (a notoriously difficult problem to solve), even a smaller system can be helpful.
> 
> Prior art: Python-MIP includes an implementation:
> https://github.com/coin-or/python-mip/blob/6044fc8f0414d71430e94b8a08573b695dc35b5a/mip/conflict.py#L15 
> 
> \### Dual feasibility report
> 
> We have \`primal\_feasibility\_report\`, but we need a way to verify the feasibility of dual solutions
> 
> This could also include a way to check whether the primal and dual objectives match, and whether \`objective\_value(model)\` matches \`value(objective\_function(model))\`. See https://github.com/jump-dev/HiGHS.jl/issues/223

> <https://github.com/jump-dev/JuMP.jl/issues/3791>
>
> I just watched @blnicho's JuMP-dev talk https://youtu.be/G1tW68vrOBM. (Thanks Be…thany! It was a very useful talk and exactly what I was after.)
> 
> It has come up in various discussions over the years, but one of the biggest differences between JuMP and Pyomo is how we handle namespacing.
> 
> \## JuMP
> 
> JuMP uses a single global namespace, and all objects are stored in the \`object\_dictionary(model)\`:
> \`\`\`julia
> model = Model()
> @variable(model, x)
> model\[:x\] === x # true
> @variable(model, x) # Errors. x already exists.
> \`\`\`
> The downside is that for unrelated model parts to compose, you need to use anonymous variables, but this has poorer printing, and the macro syntax is more limited.
> \`\`\`julia
> model = Model()
> x = model\[:x\] = @variable(model, base\_name = "x")
> model\[:x\] === x # true
> x = model\[:x2\] = @variable(model, base\_name = "x")
> \`\`\`
> 
> \## Pyomo
> 
> Pyomo has blocks. At a high level, blocks create a new namespace, and can be nested inside each other. You can delete/deactive entire blocks as a single component.
> 
> In Julia-land, we might do:
> \`\`\`julia
> model = Model()
> @block(model, blocks\[1:2\])
> @variable(blocks\[1\], x)
> blocks\[1\]\[:x\] === x # true
> @variable(blocks\[2\], x)
> blocks\[1\]\[:x\] === x # false: Julia binding \`x\` has been replaced in this scope
> blocks\[2\]\[:x\] === x # true
> @objective(model, Min, sum(b\[:x\] for b in blocks))
> \`\`\`
> 
> \## Next steps
> 
> I don't have a concrete syntax proposal, or believe that we should necessarily implement this, but I'd be open to exploring possibilities, and I'd very much like to come up with a few examples where JuMP's current syntax is limiting and the block/namespace would be beneficial.
> 
> I don't think this would change anything at the MOI level. This is strictly a JuMP-level feature.
> 
> I also don't know if this is that easy to explore in a JuMP-extension, but maybe it is. A prototype could store everything in \`model.ext\[:block\_extension\]\` for now.
> 
> Some questions:
> 
> \* What exactly is the struct of \`block\[1\]\`?
> \* Are blocks added to the model lazily or eagerly?
> \* You could imagine a block being a self-contained JuMP model with no solver backing it, and when \`block\` is referenced in the \`model\` scope everything gets copied across.
> \* Dropping a breadcrumb to https://github.com/geoffleyland/rima and https://github.com/geoffleyland/rima/blob/master/papers/Rima-ORSNZ-paper-2010.pdf

---

<div class="post-metadata">

**Author:** ![slwu89](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/slwu89/32/217323_2.png) [@slwu89](https://discourse.julialang.org/u/slwu89)\
**Post date:** [September 4, 2024, 3:50pm UTC](https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303/7 "2024-09-04T15:50:41Z")

</div>

Thanks! Yeah, regarding the nonlinear rewrite, I just meant to ask if the resolution of that work led to any more time to think about this stuff.

Thanks for the issue links, the one regarding namespaces is quite interesting. I have a few questions:

- Has there been any examples people have come up with (or questions you have seen here) where the current syntax is limiting and the pyomo style blocks would help?
- I’m aware this is all very conceptual at the moment, but if a constraint that included variables in 2 different blocks was used, where would it be stored?
- Do you see any relationship between the namespace ideas issue and this one? [Improve support for relational algebra · Issue #3438 · jump-dev/JuMP.jl (github.com)](https://github.com/jump-dev/JuMP.jl/issues/3438)

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [September 4, 2024, 11:06pm UTC](https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303/8 "2024-09-04T23:06:11Z")

</div>

Answers:

- Nope. I’m very much open to suggestions! It might be a chicken and egg situation. JuMP models are built relatively “flat” because we have only a single common namespace. If we had nested namespaces, people might write their models differently… I think it’s more a thing of: blocks seem to work quite well for Pyomo. Perhaps they could work well for JuMP. Rather than a pressing need that people in the community are asking for.

- In a common parent block of the two variables (however high up that might be)

- Nope. That’s really about efficient ways to write summations over sparse index sets.

---

<div class="post-metadata">

**Author:** ![slwu89](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/slwu89/32/217323_2.png) [@slwu89](https://discourse.julialang.org/u/slwu89)\
**Post date:** [September 6, 2024, 2:45pm UTC](https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303/9 "2024-09-06T14:45:35Z")

</div>

Regarding blocks, I am certainly interested in namespacing variables. The AlgebraicJulia people have needed this capability (`part1.x`, `part2.x`, etc) for multi-physics modeling before. I’ll have to take some time to understand what they’ve done to see if its relevant. In the meantime I can watch the pyomo talk from JuMP con and see if I can somehow get a copy of this chapter [Structured Modeling with Blocks | SpringerLink](https://link.springer.com/chapter/10.1007/978-3-319-58821-6_8) to see how these are used in practice.

The support for relational algebra is interesting. How much of that is in scope for JuMP? I saw the notes about parsing and transforming into conjunctive normal form, that could get very complex depending on how many optimizations for speed to do. The dataframes syntax from the “ijklm” model looks really clear to me, and DataFrames is really fast already. I have some examples I’ve cooked up on using Catlab to do those manipulations that I’m trying to work into a blog post.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [September 8, 2024, 10:48pm UTC](https://discourse.julialang.org/t/hierarchical-compositional-construction-of-large-optimization-models/100303/10 "2024-09-08T22:48:37Z")

</div>

The Pyomo book is available on [osti.gov](http://osti.gov): [Pyomo - Optimization Modeling in Python 3rd Ed. (Book) | OSTI.GOV](https://www.osti.gov/servlets/purl/1771935)