HELD Implementation

[branch171]

HELD Implementation with vscode example cases. Note this current only implements tp_flash_impl in HELD_flash.jl and needs integrating into the flash code particularly tpflash and tpflash2 in PT.jl.

tp_flash_impl returns beta,xp,vp,Gsol which can be placed in the flash returns for tpflash and tpflash2

beta,xp,vp,Gsol = Clapeyron.tp_flash_impl(model,p,T,z, HELDTPFlash(verbose = verbose))

Please test and feed back any comments or improvements

[branch171]

The vscode examples are out of date and need the tp_flash_impl added I will push an update later.

[longemen3000]

Suggested change

	pure = split_model.(model)
	pure = split_pure_model(model)

[longemen3000]

Suggested change

	p = Vector{Float64}(undef,n)
	p = Vector{Base.promote_eltype(x,lb,ub)}(undef,n)

[longemen3000]

Hello,

Some comments.

1. i saw that a lot of the HELD code is for a box-constrained generic trust region solver, is it possible to move that code inside modules/Solvers.jl ? you can create your own file there. i would like it to wrap inside a generic solver accessible via Solvers.optimize

2. This is something i was talking to @pw0908 some time ago, JuMP is a quite heavy dependency to add, and while i don't feel really comfortable, i understand that this is the only way to implement HELD at the moment.

3. there are some performance optimizations missing, but nothing that should block this PR.

[branch171]

Yes I think the trust region stuff could easily go into the solvers.

But note there is some specific structures to the projections for HELD and Gibbs's that are unique to the HELD code. If those are defined in the HELD code then it will work.

As for JuMP if there is a way to change it that would be fine. We really only need a simplex method. I used a canned code from numerical recipes c++ that worked

Any advice on optimisations appreciated.

[longemen3000]

some specific optimization tips, some specific to Clapeyron, others more general:

• there is a macro (defined in Clapeyron) named @sum, that is equivalent to a for-loop (@sum(x[i]*exp(y[i])).
• for gradients and hessians, you could use DiffResults.jl, or DifferentiationInterface.jl. so, instead of calling f,g and h, one should instead call a function fgh, that calculates the value, gradient and hessian in one pass.
  The latter one is not installed, but it is a light dependency.
• x[1:nc] allocates, use @view(x[1:nc) instead, that will create a view into the parent array.

Also, the specific return type for tp_flash_impl was changed, to return a Clapeyron.FlashResult. (this is to support other flash types). Before we had other flash types, we had the 4 argument return type for tp-flash, and for compatibility reasons, we still return those. (for a future breaking release, tp_flash2 will be renamed to tp_flash).

[pw0908]

It looks like you can implement HiGHS without using JuMP:
https://ergo-code.github.io/HiGHS/dev/interfaces/julia/

Not as pretty, but more reasonable?

[branch171]

I will try to use the HiGHS API and remove JuMP
And make the optimisation recommendations

I also found an annoying issue. I tested something that I thought as simple and it failed. This was air and water at 1.5 bar and 20 deg C. After a lot of frustration it turned out HELD could only solve it as water first then the air components. If the system was set up as air first then water at a certain moderately low water fraction the solver fails. It's comes up with an invalid lambda solution from the HiGHs solver that prevents the global solution. Once the components are rearranged in increasing mole fraction order it all works as expected.

Finally the trust region solver loops is prefect for parallel solution. I'd like to implement that but it can be a future update.

[pw0908]

This is something I've noticed with our other solvers as well. When you have two phases with very low concentrations, it typically prefers solving for the compositions with the order of the species in a certain order. Aside from giving guidance in the docs, I don't quite know how to handle this best?

[longemen3000]

You could sort the model in the following way:

ix = sortperm(z)
model_sort = each_split_model(model, ix)

An equivalent way is to, instead of doing x[nc] = 1 - sum(x), is doing this in the index of the component with the biggest composition, but that requires some extra care to handle the indices.

[branch171]

Usingthe sorted model works but so far I have only one case where that needed to be done. Would be good to understand why its needed and if there are any other test cases that have a similar issue that we can test.

For that same problematic case if you increase the water molefraction eventually it all works with no sorting.

I have found that the associating components need to come first. This seems to be robust.

However we may need some scaling as hessian can become ill conditioned for systems with small mole fractions.

[longemen3000]

probably, one way check if the model sort is even necessary is to check the initial K-values, and if, for example , !(1/n<= Kmin*Kmax <= n), then, perform the sorting. where n is just a constant (100?). (an equivalent way would be abs(log(Kmin*Kmax)) < log(n))

[branch171]

If have going through some testing that the statement in the original HELD paper that the initial M set can be proven to be all that's required to initiate the method is practically not true. There are some cases where for larger component sets this does not provide a good enough cutting plane. Also in the early phase of the solution when UB and LB are far apart you need all the LB solutions below the current UB otherwise you can get very poor cutting planes that case the method to fail. Once things are tighter to the approaching hyper plane solution the stated algorithm works robustly.

[branch171]

Ok what I have found is during the OPxv phase the Lambdas are very loose early on. As such for small mole fraction components you can get a sudden large Lambda that prevents the IPxv phase from finding a minimum.

The solution is to provide an upper and lower bound on the lambdas in the OPxv phase. This works but I have no way yet of knowing how to set this at the start. Certainly (G0 -Gi)/(x0 - xi) for the initial M set is a poor estimate. It's easy once I know the answer.

[longemen3000]

Hello,

I saw in some of your recent commits that you had some problems with the density solver used by GERG2008, is that right? i was just debugging #376, and there was indeed a problem with the initial value suggestions, but I'm interested if the cases are still failing (to add those cases into the test infrastructure)

[branch171]

Hi

I stopped using the Clapeyron volume solution and made my own.

The issue I had was water at 28 barg and 30 deg C. The Clapeyron solver could not identify the correct volume.

It turns out you need to identify the spinodal points and use these as criteria to select the correct root.

The issue is the many Maxwell loops with water and selecting what looks like the minimum Gibbs and mechanically stable solution but they are inside the spinodal metastable region.

I will post the test.

[branch171]

Also I have made many improvements to HELD Code since my original PR. Not sure if you get these updates. Let me know if I need to do another PR so you can get them.

The main one is moving the Trust Region to Clapeyron solvers and adding a trust region sub-problem solver that is more robust. It was make this update that I saw the issue with GERG2008 volumes.

[branch171]

This is code that fails to solve

using Clapeyron

components = ["carbon dioxide","nitrogen","water"]
model = GERG2008(components)

p = 64.0e5
T = 30+273.15

zs=[0.4975080785711593, 0.0049838428576813995]
z=append!(zs,1.0 - zs[1] - zs[2])

vol = volume(model,p,T,z)
mw = Clapeyron.molecular_weight(model,z)
println("mw = $(mw)")
rho = mw/vol
println("Volume = $(vol) rho = $(rho)")
rhom = mass_density(model,p,T,z)
println("rhom = $(rhom)")

gives

mw = 0.030997443505108243
Volume = NaN rho = NaN
rhom = NaN

[branch171]

The answer should be

Volume by HELD = 3.710503849983332e-5 m3/mol
rho by HELD = 835.3971524715569 kg/m3

[longemen3000]

It was absolutely related, it works with with the most recent commit. an easier way to test this is to check if the value of Clapeyron.x0_volume(model,p,T,z,phase = :liquid) is lower than volume(model,p,T,z,phase = :l). if an initial guess for the liquid volume is outside the permissible range, empiric models tend to fail. (EoS with more theoretical bases tend to have an spinodal that allows better convergence properties)

Before:

julia> Clapeyron.x0_volume(model,p,T,z,phase = :l)
4.614421464979334e-5

After

julia> Clapeyron.x0_volume(model,p,T,z,phase = :l)
2.3847375780049338e-5

[branch171]

Just as an observation for GERG2008. If you try to look at water at say 30 degrees C and 28 bar you will find the most stable volume including mechanical stability is a false root, at least from the reality of the required solution. This is because for compressed liquids the spinodal vapour root eventually disappears. For GERG2008 they probably don't care as it's a natural gas eos. But it's a pain if you want consistency across all pure components in a mixture. It means you can't dehydrate CO2 with compression and cooling. BTW false volume solutions cause a major problem to held as the initial M set is not on the true surface so the hyperplane solution comes up with a bad lambda solution. Saft eos,s seem to work well. My only issue is they are poor near the critical region and for CO2 this is important as we need to correctly predict near critical and super critical phase behaviour. I'm hoping for a cross over correction to SAFTVRMie may become available.

[branch171]

Note the Clapeyron Volume solution seems to be OK. So its my solution method that is getting the wrong root.

[longemen3000]

hmmm, if i remember correctly, false roots are vapour phases, right? if that is the case, and you are far from the mechanical critical point, maybe this criteria could help:

function return_true_root(model,v_liquid,v_vapour,T,x)
  gl = Clapeyron.VT0.gibbs_free_energy(model,v_liquid,T,x)
  gv = Clapeyron.VT0.gibbs_free_energy(model,v_vapour,T,x)
  if gv < gl
    B = second_virial_coefficient(model,T,x)
    if v_vapour < -2B #approximation of the spinodal, point when second order virial model diverges
      #for low temps (< 0.9Tc) this should be true. for high temps, false roots tend to disappear.
      return v_liquid
    else
      return v_vapour
    end 
  else
    return v_liquid
  end
end

Another way would be to approximate the gas gibbs energy by a virial model, check if the gibbs energy of the tentative phase is within reasonable bounds.

[pw0908]

For the simplex method, I found this in julia: https://github.com/davidagold/simplex.jl/blob/master/simplex.jl

Looks abandoned and coded in a very interesting way. Could we adapt it?