Feasibility Inquiry: Building a Chemical Process Simulation and Optimization Package with JuMP

Specific Domains
,
1

I am a user of JuMP and am aware that the Institute for the Design of Advanced Energy Systems (IDAES) is a recently developed tool built on Python that utilizes the Pyomo optimization modeling language (https://doi.org/10.1002/amp2.10095) for process simulation and optimization. I find IDAES leverages the block-oriented modeling capability of pyomo. I want to know the possibility of building a similar software based on JuMP. Note that this is a vast and intricate engineering project, and I am particularly interested in understanding the possible challenges and advantages of using JuMP. Such as performance and complexity issues.

1 Like
2

This is a fairly broad topic. @pulsipher might have some insight.

Some high level thoughts:

  • yes, you could build something similar in JuMP
  • with care, performance could be similar or better

The main consideration is that JuMP doesnt have Pyomos block structure, so you shouldnā€™t look to replicate the design of the code. Youā€™d need to rethink how to build it in a ā€œJulianā€ fashion.

Some projects to look at

  • InfiniteOpt.jl
  • PowerSimulations.jl

See also

2 Likes
3

We have an ongoing project in the Julia Lab (+ collaborations elsewhere) to build a chemical process simulation library in Julia. I am hoping it will be revealed at a football stadium this summer :wink: , but weā€™ll see about that.

Itā€™s being developed using ModelingToolkit instead for a variety of reasons.

  1. Normally the best way to solve a DAE is not to solve the equations that the user has specified. Thatā€™s a long and deep topic with lots of little algorithms that all add up, and itā€™s the bread and butter of what ModelingToolkit does.
  2. ModelingToolkit supports acausal component-based modeling, which is a generalization of the ā€œblock-oriented modelingā€ and allows for a lot of symbolic simplification. The reasons and design is similar in spirit to something like Modelica, and so this work of a chemical process simulator on OpenModelica is a good background as to why this would be used in comparison to say Aspen+ (https://pubs.acs.org/doi/epdf/10.1021/acs.iecr.9b00104).
  3. The best way to do the forward simulation of such equations is generally not to treat it as an optimization problem, rather to simulate the DAEs or directly solve the nonlinear system, depending on whether itā€™s a dynamic or steady state problem. So youā€™d generally want to make use of DifferentialEquations.jl or NonlinearSolve.jl in most cases (or BoundaryValueDiffEq.jl).
  4. A fully implicit discretization in time does not have error control over the truncation error, which in term means that the solution has a non-trivial error with respect to the chosen dt that is not estimated or shared with the user.
  5. For inverse problems where you are doing an optimization, itā€™s problem-dependent whether doing the discretization in time as part of the optimization is good or bad. There are examples that go both ways, and so having the flexibility to choose to treat it via DAE solvers or fully implicit can help (or the mixture, i.e. as a multi-point BVP).
  6. Importantly, ModelingToolkit comes with a standard library of other components so then you can mix and match the chemical process simulation with other aspects, like a hydraulic system or a power simulation. The core of an acausal modeling system is to help bring different models together.

There are still some downsides, which are being worked out as part of the project. For one thing, we need simpler connections to multi-point BVPs, so while BoundaryValueDiffEq.jl just got a revamp over the last year, ModelingToolkit doesnā€™t target it effectively. Also, if nonlinear constraints are added then it does need to fall back to a form of feasibility problem that uses optimization, and the codegen to optimizers is currently not as optimized as JuMP (though thereā€™s ongoing work on Optimization.jl and targeting of JuMP code directly from MTKā€™s symbolic form).

2 Likes
4

So, do you mean the software IDAES (bulid up on python and pyomo) is not a good choice for chemical process simulation and optimization, using Modelica is better?

5

I think thereā€™s good reasons to take an acausal modeling approach to the development of the system. The OpenModelica results, while not ending up in an actively maintained repository, do show some of the benefits of doing so. And gPROMs is effectively an acausal modeling system as well and I think is a nice and clear demonstration of how acausal modeling can help chemical process simulation be more robust.

Of course, there are trade-offs. The main trade-off being that an acausal modeling system can be wildly more complicated to engineer: simply using the equations given to you by the user is the most straightforward thing to do. Also, it can take a lot of design to get the connection to causal systems right, which is required for building control systems. The other is that the optimization-based approach more naturally lends itself to feasibility solves that can enforce nonlinear constraints. While itā€™s feasible in a dynamical approach (differential inclusions), itā€™s definitely not straightforward, and thus any handling of these cases would need to fall back to optimization tooling in some form (or multi-point BVP). But you can still apply the equation transformations to an extent so it can still be useful to mix these techniques with an optimization approach.

2 Likes
6

Thank you for your enthusiastic reply. Iā€™ll continue to explore various advanced modeling frameworks. It seems that IDAES is also based on an acausal modeling approach (perhaps I misunderstood your statement; I consider the equation-oriented modeling method as the acausal modeling approach).

7

Having closely worked with the IDAES team, the main benefit of developing the process simulation framework in an AML environment like Pyomo is being able to conduct advanced analysis and design techniques (e.g., robust optimization, flexibility analysis, stochastic optimization, etc.) which require constrained programming. For these types of problems, the team and myself found that shooting methods (i.e., forward simulating the DAEs) didnā€™t perform well against simultaneous methods (e.g., discretize then optimize).

Depending on your goals, it would be possible in principle to make something analogous to IDAES in Julia. Perhaps using the SciML framework if simulation and unconstrained optimization is your primary motivation, or using JuMP if you have similar objectives to the IDAES team. However, it would be an enormous undertaking. It has taken 40+ IDAES team members 6+ years to get it to where it is now. As a ChemE researcher and a Julia fanboy, I would love something like this to exist but it seems like too large an undertaking, especially considering that the IDAES framework already exists.

I will also note that one fundamental challenge the IDAES team has is that good open-source thermodynamic property data is few and far between. This definitely gives commercial tools like Aspen a big advantage.

3 Likes
8

Thanks, sincerely.

2 Likes
9

Hey Chris!

You cheekily mentioned a big reveal at JuliaCon2024ā€¦ its soon July and Iā€™m deeply invested in the status of this library!

My questions:

  1. Is it on track to be revealed (and released) at JuliaCon?
  2. If it is delayed, where can I follow and or contribute to its development?

Iā€™m currently implementing process unit models in pyomo and julia for opensource, and a library like this would be very useful indeed.

We will make sure to give the library a shoutout to our followers at PolyModels Hub!
https://www.linkedin.com/company/polymodels-hub/?originalSubdomain=uk

Cheers,
Sam

1 Like
10

Itā€™s both going to be shown but also a bit delayed.

Whatā€™ll be shown is an early flow sheeting system. No GUI integration yet, also missing some optimizations, but the bones are there.

1 Like
11

It feels like my prayers have been answeredā€¦ it even interfaces with Clapeyron.jl !

Thank you for the great work, I canā€™t wait to build on top of this. It will certainly unlock a lot of value and makes it possible for the community to build a large library of process unit models, which is sorely lacking.

12

How should I participate in the development of this project now? I have nearly 5 years of experience in process simulation software development. I am familiar with sequential unit equipment algorithms and EO method modeling. I look forward to contributing.

13

Itā€™s not open source yet, but it should be soon?

14

Thank you for your interest in this initiative. Weā€™re currently focused on fixing issues to have enough features for public release. One of the key challenges weā€™re addressing is the integration of rigorous thermophysical property calculations within the governing unit operation equations within ModelingToolkit.

Given your long experience in process simulation, would you have insights into how existing solutions handle this integration? Your perspective could be help in fixing this part and accelerate the public release.

15

According to the experience of gPROMS and our own development experience, the thermodynamics library should be implicitly called as an external module of the model equation, such as this call statement Y(,j)*phys_prop.VapourFugacityCoefficient(T(j),P(j),Y(,j)) = X(,j)*phys_prop.LiquidFugacityCoefficient(T(j),P(j),X(,j));. During development, the return value and derivative value interface of the thermodynamics module should be developed, and then registered using the external function function of mtk. I am currently studying the integration of multiflashā€™s python interface into mtk. If there is any progress, I will be happy to share it with your team.multiflash

1 Like
16

Thanks for sharing your experience. From what you described, it seems to be close to what we are doing with clapeyron.jl developers.

But, I think it would be great to have the possibility to interface with other libraries including proprietary ones for thermophysical property. With the fixes I mentioned, we should be able to put a decent code base open for further developments like this.

1 Like
17

I would like to ask some technical questions. I am trying to connect clapeyron to mtk, but I have some systemic problems when calling it. Can you give me some suggestions for connection?

18

How are you performing the registration?

19

If you have any questions about Clapeyron, feel free to reach out for help, or leave an issue in the repository

1 Like
20

I dont know what it the right way to registration, I used this code

@variables  T x[1:numcomp] y[1:numcomp]
@parameters  Tc[1:numcomp] Pc[1:numcomp] Ļ‰[1:numcomp]  z[1:numcomp] vaporFraction P
function simpleK(i, T, P, Tc, Pc, Ļ‰)
    return Pc[i] * exp(5.37 * (1 + Ļ‰[i]) * (1 - Tc[i] / T)) / P
end
function CalcK(T,P,x,y)
    K = Clapeyron.fugacity_coefficient(model,P,T,x,phase=:l)./Clapeyron.fugacity_coefficient(model,P,T,y,phase=:v)
    return K
end
function LiquidFug(T,P,x)
    return Clapeyron.fugacity_coefficient(model,P,T,x,phase=:l)
end
function VaporFug(T,P,x)
    return Clapeyron.fugacity_coefficient(model,P,T,x,phase=:v)
end
function simpleK(i, T, P, Tc, Pc, Ļ‰)
    return Pc[i] * exp(5.37 * (1 + Ļ‰[i]) * (1 - Tc[i] / T)) / P
end
function CalcK(T,P,x,y)
    K = Clapeyron.fugacity_coefficient(model,P,T,x,phase=:l)./Clapeyron.fugacity_coefficient(model,P,T,y,phase=:v)
    return K
end
function LiquidFug(T,P,x)
    return Clapeyron.fugacity_coefficient(model,P,T,x,phase=:l)
end
function VaporFug(T,P,x)
    return Clapeyron.fugacity_coefficient(model,P,T,x,phase=:v)
end
eqs = [
    sum(y) ~ 1,
    [z[i] ~ x[i]*(1-vaporFraction)+y[i]*vaporFraction for i in 1:numcomp]...,
    [y.~ x.* VaporFug(T, P,x)]...,
]
@mtkbuild ns = NonlinearSystem(eqs, [T, x...,y...,], [Tc..., Pc..., Ļ‰...,z...,vaporFraction,P])

I got the error about

 #= C:\Users\.julia\packages\Symbolics\2UpZj\src\arrays.jl:591 =#
substitute(j, dict)end), ~(~i))...,)) on expression (broadcast(*, x, broadcast(/, broadcast(exp, broadcast(/, broadcast(/, Clapeyron.VT_partial_property(RK{BasicIdeal, SoaveAlpha, NoTranslation, vdW1fRule}("ethane", "ethylene", "propylene", "methane", "propane"), Clapeyron._volume(RK{BasicIdeal, SoaveAlpha, NoTranslation, vdW1fRule}("ethane", "ethylene", "propylene", "methane", "propane"), P, T, x, v, true, nothing), T, x, eos_res), 8.31446261815324), Ref(T))), Ref((-0.12027235504272604Clapeyron.āˆ‚fāˆ‚V(RK{BasicIdeal, SoaveAlpha, NoTranslation, vdW1fRule}("ethane", "ethylene", "propylene", "methane", "propane"), Clapeyron._volume(RK{BasicIdeal, SoaveAlpha, NoTranslation, vdW1fRule}("ethane", "ethylene", "propylene", "methane", "propane"), P, T, x, v, true, nothing), T, x)*Clapeyron._volume(RK{BasicIdeal, SoaveAlpha, NoTranslation, vdW1fRule}("ethane", "ethylene", "propylene", "methane", "propane"), P, T, x, v, true, nothing)) / (T*Symbolics._mapreduce(identity, +, x, Colon(), (:init => false,)))))))[i]

Stacktrace:

Can you give me some examples?"