Michelsen Multiphase preliminary work

[lucpaoli]

Hi, this is my initial work on a multiphase michelsen algorithm. The only real change is the Rachford-Rice equation, which is reformulated to a convex minimisation problem, as in

(1) Michelsen, M. L. Calculation of Multiphase Equilibrium. Computers & Chemical Engineering 1994, 18 (7), 545–550. https://doi.org/10.1016/0098-1354(93)E0017-4.

However, rather than following the algorithm detailed by him, this instead relies on interior point newton as a local optimiser - defined within Optim.jl. I have previously worked with a manual implementation, but the control structure for this is annoying and it does not perform as well as the interior point newton within Optim. I'm not sure if the interior point newton is defined within NLSolvers.jl, which I think you try to keep as your sole dependency for optimisation?

The multiphase RR correctly reproduces Tables 1 and 3 within the paper above, though Table 2a and 2b do not work for reasons I cannot tell.

The actual algorithm beyond the multiphase RR is much the same, though it operates on the mole number matrix, x, rather than a set of K-values. This is because as far as I know, when working with K-values, care has to be taken that the selected reference phase does not disappear during the multiphase RR calculation.

The next step from this in my view is to get stability analysis working properly on the resulting phase distribution, as that can then be used to add phases and recalculate the flash. The primary difficulty here is avoiding convergence to a trivial solution. I was thinking it could also be interesting to combine this with the HELD stability analysis.

I attempted to use your tpd() function, however this appears to return incorrect minima, which result in converging to the wrong phase distribution. This was tested using the system

model = EPPR78(["methane", "hydrogen sulfide"])
p = 4.052e6
T = 180.0
z = [0.5, 0.5]

and verified that the value returned by tpd() is not a minima of the full TPD function, or the equivalent tm function. I'm not certain what the issue is - perhaps that's just initial guesses, though the extreme values return suggest perhaps not?

I'm not creating this PR to immediately merge, but more to make you aware of the implementation and invite any comments / modifications that would make this more production-ready. I cannot test the implementation Michelsentp_flash.jl, as this appears to error completely on my machine, but I do not believe this allows for multiphase calculations.

[codecov-commenter]

Codecov Report

Base: 0.00% // Head: 0.00% // No change to project coverage :thumbsup:

Coverage data is based on head (5458e2e) compared to base (ebb810f). Patch coverage: 0.00% of modified lines in pull request are covered.

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #108    +/-   ##
=======================================
  Coverage    0.00%   0.00%            
=======================================
  Files         209     211     +2     
  Lines       13632   13897   +265     
=======================================
- Misses      13632   13897   +265     

Impacted Files	Coverage Δ
src/methods/methods.jl	0.00% <0.00%> (ø)
src/methods/pT.jl	0.00% <0.00%> (ø)
...ethods/property_solvers/multicomponent/tp_flash.jl	0.00% <0.00%> (ø)
...multicomponent/tp_flash/GuptaMultiphasetp_flash.jl	0.00% <0.00%> (ø)
...icomponent/tp_flash/MichelsenMultiphasetp_flash.jl	0.00% <0.00%> (ø)
...rty_solvers/singlecomponent/saturation/ChemPotV.jl	0.00% <0.00%> (ø)
src/methods/property_solvers/volume.jl	0.00% <0.00%> (ø)
src/methods/stability.jl	0.00% <0.00%> (ø)
src/methods/tpd.jl	0.00% <0.00%> (ø)
src/solvers/Solvers.jl	0.00% <ø> (ø)
... and 4 more

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

:umbrella: View full report at Codecov.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.

[longemen3000]

I was looking at the IPNewton use, and it's mainly for Newton with bounds, right? (in the 0x1 hypercube). in that case, we can use optimize(f,x0,NLSolvers.ActiveBox(;bounds)) (or something like that) that performs a projected gradient method. in theory they should perform the same.

I cannot test the implementation Michelsentp_flash.jl, as this appears to error completely on my machine, but I do not believe this allows for multiphase calculations.

if i remember correctly, there there was a bug, that was fixed on https://github.com/ypaul21/Clapeyron.jl/pull/98/commits/9885fbf578d839b40af7170a85f9f8568891f891 . and yes, it is only for 2-phase calculations.

On tpd, i have no idea why it is failing.

[longemen3000]

on moving all solvers to NLSolvers, i did some limited testing: Optim.IPNewton vs NLSolvers.ActiveBox (with positive cholesky factorization)

* Status: success

 * Candidate solution
    Final objective value:     -1.395441e+00

 * Found with
    Algorithm:     Interior Point Newton

 * Convergence measures
    |x - x'|               = 4.70e-10 ≰ 0.0e+00
    |x - x'|/|x'|          = 9.40e-10 ≰ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 8.09e-11 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    16
    f(x) calls:    43
    ∇f(x) calls:   43

Results of minimization

* Algorithm:
  ActiveBox

* Candidate solution:
  Final objective value:    -1.40e+00
  Final gradient norm:      1.51e-14
  Final projected gradient norm:  1.51e-14

  Initial objective value:  -1.40e+00
  Initial gradient norm:    3.09e-04

* Convergence measures
  |x - x'|              = 8.70e-08 <= 0.00e+00 (false)
  |x - x'|/|x|          = 1.23e-07 <= 0.00e+00 (false)
  |f(x) - f(x')|        = 7.99e-15 <= 0.00e+00 (false)
  |f(x) - f(x')|/|f(x)| = 5.73e-15 <= 0.00e+00 (false)
  |x - P(x - g(x))|     = 1.51e-14 <= 1.00e-08 (true)
  |g(x)|                = 1.51e-14 <= 1.00e-08 (true)
  |g(x)|/|g(x₀)|        = 4.89e-11 <= 0.00e+00 (false)

* Work counters
  Seconds run:   6.40e-01
  Iterations:    2

NLSolvers.ActiveBox seems to be faster in that regard. it was expected, as Optim.IPNewton also performs other things to be able to accomodate non-linear constraints.

NLSolvers.Linesearch(Newton(),Dogleg) vs NLsolve(method = :trust_region) results seems to be the same.

NLSolvers.Anderson vs NLsolve(method = :anderson) i did have a problem with the NLSolvers.jl package https://github.com/JuliaNLSolvers/NLSolvers.jl/issues/44 . also i got bit by https://github.com/JuliaNLSolvers/NLSolvers.jl/issues/25 . the latter can be fixed on our part, the former requires upstream support, so for the moment, NLsolve it is

[lucpaoli]

Hi Andreas, yup the only reason to use the IP Newton is for a bounded optimiser. The bounds imposed are just linear, for the purpose of keeping betas in physical range, [0, 1]. As the problem is convex, convergence should not be an issue beyond that. For some reason, the hand coded Newton solver had some issues - for example, running into a singular Hessian on occasion, or otherwise getting "stuck".

The projected gradient method looks great, sorry for not working on this more.

[gustavochm]

Hi both!

I recently attended to the Advanced thermodynamic course and during the course I coded this multiphase flash for Clapyeron. I actually got a running version of the phase split calculation, then phase equilibria by ASS and finally Gibbs energy minimization if necessary.

For now, the code is pretty "ugly" but the good thing is that handles activation/elimination of phases on both ASS and the energy minimization.

@lucpaoli, if you are still into this we could use that to code a "nicer" working version to Clayperon :)

Gustavo

[pkofod]

https://github.com/JuliaRegistries/General/pull/74366

[longemen3000]

@lucpaoli there was a bug on TPD, now seems to give more sensical phases:

tpd(model,p,T,z)
([[0.059765812629097095, 0.9402341873709029], [0.9396783613334813, 0.060321638666518726]], [-0.15088179030413396, -0.14903526879769657], [:liquid, :liquid], [:vapour, :vapour])

[pw0908]

@longemen3000 What's left to do here?

[longemen3000]

I would rebase this entire branch, using the Anderson solver in cDFT

[longemen3000]

replaced by #272