search

ClapeyronThermo / Clapeyron.jl

Clapeyron provides a framework for the development and use of fluid-thermodynamic models, including SAFT, cubic, activity, multi-parameter, and COSMO-SAC.
MIT License
188 stars 47 forks source link

Parameter Estimation for sle_solubility #224

Open samuel-zhang01 opened 8 months ago

samuel-zhang01 commented 8 months ago

Amazing devs of Clapeyron.jl,

Sorry to trouble you guys again. I'm exploring parameter estimation for PCSAFT of APIs. Due to the low vapour pressure of APIs, obtaining $p{sat}$ and $\rho{l,sat}$ experimentally is quite challanging. experimental sle_solubility data of APIs in different solutes though should be a good substitute to do parameterisation. Is there a built-in support for param estimation for segment, sigma, epsilon, epsilon_assoc, and bondvol?

# Here is an imaginary param estimation code that uses sle_solubility
# init
using Pkg
Pkg.activate("..")
using Clapeyron, BlackBoxOptim

# Defining a sle_solubility function
function EXPsolubility(model::CompositeModel,T)
    sol = sle_solubility(model,101325.,T,[1.,0.];solute=[solute1])
    return sol[2]
end

# Defining solute and solvent, will be fetched from the database later
solute1="nicotinamide2016"
solvent1="ethanol2016"

# Defining a composite model
model = CompositeModel([solvent1,solute1];fluid=PCSAFT,solid=SolidHfus,gas=PCSAFT,gas_userlocations=["data/sle_specific"],fluid_userlocations=["data/sle_specific"],solid_userlocations=["data/sle_specific"])

# Creating the Dict for estimation
toestimate = [
    Dict(
        :param => :epsilon_assoc,
        :lower => 1000.,
        :upper => 3000.,
        :guess => 2500.
    ),
    Dict(
        :param => :bondvol,
        :lower => 0.02,
        :upper => 0.04,
        :guess => 0.03
    )
];

# Estimation
estimator,objective,initial,upper,lower = Estimation(model,toestimate,["data/EXPsolubility.csv"]);

# BBOptim
nparams = length(initial)
bounds  = [(lower[i],upper[i]) for i in 1:nparams]
bounds
result = BlackBoxOptim.bboptimize(objective; 
        SearchRange = bounds, 
        NumDimensions = nparams,
        MaxSteps=200,
        PopulationSize = 1000,
        TraceMode=:verbose)

params = BlackBoxOptim.best_candidate(result);

here is the error code:

MethodError: no method matching sle_solubility(::CompositeModel{PCSAFT{BasicIdeal, Float64}, PCSAFT{BasicIdeal, Float64}, SolidHfus, PCSAFT{BasicIdeal, Float64}, Nothing}, ::Float64)

Closest candidates are:
  sle_solubility(::CompositeModel, ::Any, !Matched::Any, !Matched::Any; solute)
   @ Clapeyron C:\Users\USER\.julia\packages\Clapeyron\Gx6LP\src\methods\property_solvers\multicomponent\solids\sle_solubility.jl:9

Stacktrace:
  [1] _broadcast_getindex_evalf
    @ .\broadcast.jl:683 [inlined]
  [2] _broadcast_getindex
    @ .\broadcast.jl:656 [inlined]
  [3] getindex
    @ .\broadcast.jl:610 [inlined]
  [4] copy
    @ .\broadcast.jl:912 [inlined]
  [5] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(sle_solubility), Tuple{Base.RefValue{CompositeModel{PCSAFT{BasicIdeal, Float64}, PCSAFT{BasicIdeal, Float64}, SolidHfus, PCSAFT{BasicIdeal, Float64}, Nothing}}, Vector{Union{Missing, Float64}}}})
    @ Base.Broadcast .\broadcast.jl:873
  [6] objective_function(estimation::Estimation{CompositeModel{PCSAFT{BasicIdeal, Float64}, PCSAFT{BasicIdeal, Float64}, SolidHfus, PCSAFT{BasicIdeal, Float64}, Nothing}}, guesses::Vector{Float64})
    @ Clapeyron C:\Users\USER\.julia\packages\Clapeyron\Gx6LP\src\estimation\estimation.jl:333
  [7] objective
    @ C:\Users\USER\.julia\packages\Clapeyron\Gx6LP\src\estimation\estimation.jl:137 [inlined]
  [8] fitness(x::Vector{Float64}, p::FunctionBasedProblem{Clapeyron.var"#objective#1115"{Estimation{CompositeModel{PCSAFT{BasicIdeal, Float64}, PCSAFT{BasicIdeal, Float64}, SolidHfus, PCSAFT{BasicIdeal, Float64}, Nothing}}}, ScalarFitnessScheme{true}, ContinuousRectSearchSpace, Nothing})
    @ BlackBoxOptim C:\Users\USER\.julia\packages\BlackBoxOptim\I3lfp\src\problem.jl:61
  [9] setup_problem(func::Function, parameters::ParamsDictChain)
...
    @ BlackBoxOptim C:\Users\USER\.julia\packages\BlackBoxOptim\I3lfp\src\bboptimize.jl:92
 [13] bboptimize
    @ C:\Users\USER\.julia\packages\BlackBoxOptim\I3lfp\src\bboptimize.jl:91 [inlined]
 [14] top-level scope
    @ c:\Users\USER\Documents\Documents\GitHub Repos\PC-SAFT-Thermo-Properties-Tool\example-calculations\Parameter-Estimation\parameter_estimation test.ipynb:4

I've linked the databases used for my calculation.

Thank you so much for your help!

data.zip

longemen3000 commented 8 months ago

Seems a bug on our part, let me see what i can do

samuel-zhang01 commented 8 months ago

@longemen3000 you are a hero, thanks!

pw0908 commented 8 months ago

Hi Samuel, I took a look at your zip file. It looks like you labelled the method incorrectly in your EXPsolubility.csv file. The method used should be the name of the function you defined above:

Clapeyron Estimator,
[method = EXPsolubility],
T,out_x

This should fix the issue.

pw0908 commented 8 months ago

To address your question regarding estimating the other parameters, you can take a look at the example notebook. Since you are estimating parameters for a mixture, you'll need to specify the indices. Here is how I would code it up:

using Clapeyron, BlackBoxOptim

function EXPsolubility(model::CompositeModel,T)
    sol = sle_solubility(model,101325.,T,[1.,0.];solute=[solute1])
    return sol[2]
end

# Defining solute and solvent, will be fetched from the database later
solute1="nicotinamide2016"
solvent1="ethanol2016"

# Defining a composite model
model = CompositeModel([solvent1,solute1];fluid=PCSAFT,solid=SolidHfus,gas=PCSAFT,gas_userlocations=["data/sle_specific"],fluid_userlocations=["data/sle_specific"],solid_userlocations=["data/sle_specific"])

# Creating the Dict for estimation
toestimate = [
    Dict(
        :param => segment,
        :indices => 2,
        :lower => 1.,
        :upper => 5.,
        :guess => 2.5
    ),
    Dict(
        :param => :sigma,
        :indices => (2,2),
        :recombine => true,
        :lower => 2.8,
        :upper => 3.8,
        :guess => 3.3,
        :factor => 1e-10
    ),
    Dict(
        :param => epsilon,
        :indices => (2,2),
        :recombine => true,
        :lower => 250.,
        :upper => 500.,
        :guess => 350.
    ),
    Dict(
        :param => :epsilon_assoc,
        :indices => 4,
        :recombine => true,
        :lower => 1000.,
        :upper => 3000.,
        :guess => 2500.
    ),
    Dict(
        :param => :bondvol,
        :indices => 4,
        :recombine => true,
        :lower => 0.02,
        :upper => 0.04,
        :guess => 0.03
    )
];

# Estimation
estimator,objective,initial,upper,lower = Estimation(model,toestimate,["data/EXPsolubility.csv"]);

# BBOptim
nparams = length(initial)
bounds  = [(lower[i],upper[i]) for i in 1:nparams]
bounds
result = BlackBoxOptim.bboptimize(objective; 
        SearchRange = bounds, 
        NumDimensions = nparams,
        MaxSteps=200,
        PopulationSize = 1000,
        TraceMode=:verbose)

params = BlackBoxOptim.best_candidate(result);

I will quickly point out that fitting 5 parameters using 5 data points is very ill advised. I would recommend increasing the size of the dataset (wider temperature range or more solvents) or setting some initial values for the parameters and adjust a smaller subset of parameters (typically, people just fit segment, sigma and epsilon, setting epsilon_assoc and bondvol to be some constant).

pw0908 commented 8 months ago

@longemen3000 I realised while doing an initial run of this that, in the case where recombine=true, we should also recombine the cross association parameters. If Samuel were to just estimate the self-association parameters above, it wouldn't update the cross association parameter between the solvent and solute.

longemen3000 commented 8 months ago

Hmmm, that requires a patch then

samuel-zhang01 commented 8 months ago

@pw0908 Thanks for the very swift reply. You are right about the crossaccos values. Following the advice on param estimation, is there a way to set up the BBOptim to run for different solvents without having to rebuild toestimate each time? If not, how do I carry over the param estimation of one solvent to another? Thanks!

pw0908 commented 8 months ago

Hi @samuel-zhang01,

This will be a bit more challenging but should be doable within Clapeyron. The first thing youll need to do is define a model with the API you're trying to fit and the various solvents you intend to include:

model = CompositeModel(["solute","solvent1","solvent2",...]) # inlcude all the messy parameters here.

It is important to make sure the solute is the first component as this will mean that all the solute parameters will be index 1. To assemble your data, you will probably have separate spreadsheets for each solvent. Within Clapeyron, you can define which components are involved in which data set within the spreadsheet headers:

Clapeyron Estimator,,,
[method= EXPsolubility, species=solute solvent1],,,

This means that Clapeyron will automatically reduce your model within the estimation to a binary mixture where the EXPsolubility function you wrote should still work (where the composition is defined as [0.,1.] instead). You can then include all your datasets within the Estimation function.

To keep your life simple, I recommend you don't fit the association parameters (just set them to some reasonable constant). We haven't fixed the cross association problem yet so I don't think it'll be wise to do so either way.

Since I don't have the experimental data for this system, I can't write you a working example. But I do believe this should work.

samuel-zhang01 commented 8 months ago

Thank you so much @pw0908, I will give it a try

samuel-zhang-pfizer commented 8 months ago

Dear Devs of Clapeyron.jl,

I might need some more help regarding parameter estimation and the toestimate() function.

# using Pkg,Clapeyron
# Pkg.activate("..")
using Clapeyron, BlackBoxOptim
function EXPsolubility(model::CompositeModel,T)
    sol = sle_solubility(model,101325.,T,[1.,0.];solute=[solute1])
    return sol[2]
end

# Defining solute and solvent, will be fetched from the database later
solute1="lidocaine2013"
solvent1="n-hexane2013"
solvent2="n-heptane2013"
solvent3="ethanol2013"
solvent4="acetone2013"
databasePATH="data/sle_specific"
EXPdataPATH="data/ExpData"

# Defining a composite model
model = CompositeModel([solvent1,solvent2,solvent3,solvent4,solute1];
    fluid=PCSAFT,
    solid=SolidHfus,
    gas=PCSAFT,
    gas_userlocations=[databasePATH],
    fluid_userlocations=[databasePATH],
    solid_userlocations=[databasePATH]
    )

# Creating the Dict for estimation
toestimate = [
    Dict(
        :param => :segment,
        :indices => 5,
        :lower => 1., # 1, for segment term is a good start
        :upper => 9., # 
        :guess => 5.
    ),
    Dict(
        :param => :sigma,
        :indices => (5,5),
        :recombine => false,
        :lower => 1., # start from 1
        :upper => 5., # around 5, never more than 5
        :guess => 3.,
        :factor => 1e-10
    ),
    Dict(
        :param => :epsilon,
        :indices => (5,5),
        :recombine => false,
        :lower => 100., # start at around 100, 
        :upper => 500., # up to 500-600
        :guess => 155.  # around half way of lower bond, 150
    ),
];

estimator,objective,initial,upper,lower = Estimation(model,toestimate,[EXPdataPATH]);

nparams = length(initial)
bounds  = [(lower[i],upper[i]) for i in 1:nparams]

result = BlackBoxOptim.bboptimize(objective; 
                SearchRange = bounds, 
                NumDimensions = nparams,
                # MaxSteps=5,
                MaxTime = 3600*18,
                Method = :adaptive_de_rand_1_bin_radiuslimited,
                PopulationSize = 1000,
                TraceMode=:compact)
        params = BlackBoxOptim.best_candidate(result);

Questions:

Thanks again for helping!

Below are the experimental data used for my fittings. Parameter-Estimation.zip

pw0908 commented 8 months ago

Ill take a look. One thing I immediately notice: the paper you refer to for the lidocaine data uses electrolyte PC-SAFT along with reactive equilibrium. This is not accounted for at all in base PC-SAFT. Are you just trying to create a simplified version of the model? Without the effect of pH?

samuel-zhang01 commented 8 months ago

Ill take a look. One thing I immediately notice: the paper you refer to for the lidocaine data uses electrolyte PC-SAFT along with reactive equilibrium. This is not accounted for at all in base PC-SAFT. Are you just trying to create a simplified version of the model? Without the effect of pH?

Hi Pierre, thanks for the quick reply! Yes you are right, i'm ignoring the ePCSAFT dependency on pH and only trying to reproduce pure PCSAFT parameters

samuel-zhang-pfizer commented 7 months ago

ExpData.zip I've realised i forgot to convert solubility data from the mg/g to mol/mol, reuploading...

pw0908 commented 7 months ago

Hi Samuel, sorry for the late response. I am now back from conference and have been able to take a look at your script and made some modifications in the attached zip file: example.zip

Answering your questions directly:

  1. That combining rule within the example notebook refers to the combining rule between dispersive parameters, not associative. We have not yet updated Clapeyron to automatically update association combining rules during parameter estimation.
  2. The reason you're getting that error is because the vector you input in sle_solubility is [1,0] which tells the code that the first component is the solvent (in your case, the lidocaine when you put it as the first component). I've now fixed this in the zip file I've attached.
  3. As of yet, we haven't implemented temperature dependent combining rules. But, when you run the script I've attached, you'll find that the parameters end up being quite close to those obtained in literature, which seems to imply that the combining rule doesn't affect much?

Some brief comments about things I changed in the scripts: