AttainableRegionTutorial.jl

A dataset consists of MSMPR crystallizer operating conditions and associated measurement results. For the purpose of this tutorial data is assumed to consist of steady state crystal size distributions that all possess the same discretization.

Fields

• τ: Vector of residence times [s]

• T: Vector of temperatures [°C]

• Cf: Vector of feed concentrations [kg m⁻³]

• L: Vector of crystal sizes at which CSD data is available [m].

• n: Matrix of CSDs [m⁻⁴]; each column is a CSD.

Constructors

To generate artificial data from a model:

Dataset(residence_times, temperatures, feed_concentrations, 
        system::SystemSpecification = SystemSpecification(); 
        nBins = 101, maxSize = 2000e-6, noise_level = 0.02)

To directly supply measured data:

Dataset(residence_times, temperatures, feed_concentrations, particle_sizes, CSDs)

Contains a function describing kinetics of a crystallization mechanism and information on how many optimizable parameters occur in the function.

Fields

• f: Function describing kinetics

• Nparameters: Number of optimizable parameters in function

Constructors

FunctionalExpression(f::Function, Nparameters::Integer)

Examples

S is supersaturation ln(c/c) where c is solubility. T is temperature in °C. M is normalized suspension density.

Example 1: Growth rate definition

g = (S, T, p) -> p[1] * S^p[2] * exp(-p[3] / (T + 273.15)) 
growthrate = FunctionalExpression(g, 3) #3 parameters occur in g

Example 2: Nucleation rate definitions

b1 = (S, T, M, p) -> p[1] * S^p[2] * M^p[3] 
nucleationrate1 = FunctionalExpression(b1, 3) #3 parameters occur in b

b2 = (S, T, M, p) -> p[1] * S^p[2]
nucleationrate2 = FunctionalExpression(b2, 2) 

Note that function signature must be maintained despite M not appearing on r.h.s. of b2.

Example 3: Solubility line

cstar = (T, p) -> 3.79e-2*T^2 + 3.77e-1*T + 2.07e1
solubility = FunctionalExpression(cstar, 0) 

A system specification contains material constants, a solubility expression, crystallization kinetics and parameters involved in those expressions.

Fields

• density: Density of the crystalline material [kg m⁻³]

• shape_factor: Volume shape factor k, so that volume of a crystal is V = kL³

• solubility: Solubility function [kg m⁻³]

• growth_rate: Growth rate [m s⁻¹]

• nucleation_rate: Nucleation rate [m⁻³ s⁻¹]

• parameters: Vector of parameters occuring in solubility function and growth and nucleation rates

Constructors

b = (S, T, M, p) -> p[1] * S^p[2] * M^p[3]
g = (S, T, p) -> p[1] * S^p[2] * exp(- p[3] / 8.31441 / (T+273.15))
cstar = (T,p) -> 3.79e-2*T^2 + 3.77e-1*T + 2.07e1
SystemSpecification(; density = 1200.0,
         shape_factor = pi / 6,
         solubility = FunctionalExpression(cstar, 0),
         nucleation_rate = FunctionalExpression(b, 3),
         growth_rate = FunctionalExpression(g, 3),
         parameters = [3.34e-4, 1.1, 1.44e4, 3e5, 2.0, 1.6])

ps, pg, pn = __indexp(system::SystemSpecification)

Internal implementation detail. Convenience function that returns the optimizable parameters in the solubility, the growth rate and the nucleation rate (ps, pg and pn, respectively) occurring in system.

f = error_fun(parameters, (system::SystemSpecification, dataset::Dataset))

Objective function that is minimized during parameter estimation. The implementation here minimizes the sum of squared differences between simulated and experimental CSDs, n: $f = \sum_i^p\sum_j^q \left(\ln(n_{\mathrm{sim},i,j}+1) - \ln(n_{\mathrm{exp},i,j}+1)\right)^2$ where p is the number of datasets and q is the number of particle size bins.

estimatedSystem::SystemSpecification = estimate_parameters(system::SystemSpecification, 
    dataset::Dataset; initial_guess = nothing, error_fun = error_fun, algorithm = PolyOpt())

system contains the model whose parameters should be estimated, dataset is the data that should be matched by the model, initial_guess is the initial parameter guess, error_fun is a function returning the objective function value (a metric quantifying the difference between the data and the model) and algorithm is the optimization algorithm used to minimize the objective function.

initial_guess defaults to nothing (initial guess is taken within +/- 10% of the supplied parameters in system) or can be a user-supplied guess (must match length of parameters in system).

algorithm must be one of PolyOpt(), BBO_adaptive_de_rand_1_bin_radiuslimited(), or LBFGS(),

PolyOpt() is a Poly-algorithm that combines the stochastic ADAM optimizer as built into Optimisers.jl and a gradient-based quasi-Newton method (see LBFGS() below).

BBO_adaptive_de_rand_1_bin_radiuslimited() is an evolutionary algorithm from the package BlackBoxOptim.jl (see also Optimization.jl's documentation page).

LBFGS() is a quasi-Newton (Broyden-Fletcher-Goldfarb-Shanno) algorithm from the package Optim.jl

diffMT = find_Css(Css, system, τ, T, Cf)

calculates the magma/suspension density using the overall mass balance of the MSMPR crystallizer (M) and the analytical solution of the PBE (at steady state) for the MSMPR crystallizer ($M_2$) and returns the difference between the two.

In mathematical form: $M = C_f - C_{ss}$

$M2 = k_v \rho_c \int\limits_0^\infty L^3n(L) dL = 6 k_v \rho_c B τ^4 G^3$

where $k_v$ is the shape factor, $\rho_c$ is the crystal density, $B$ is the nucleation rate, $G$ is the growth rate, $C_f$ is the feed concentration, τ is the residence time, and $C_{ss}$ is the steady state concentration.

generate_attainable_region(system::SystemSpecification; temperature_range = (5.0, 75.0), 
        feed_conc_range = nothing,
        residence_time_range = (600.0, 36000.0), npoints_temperature = 200, npoints_feed_conc = 200,
        npoints_residence_time = 200)

TODO: complete docstring.

plot_CSDs(dataset, system = nothing; ids = 1:length(dataset.T), mode = :combined, 
        figure = (fontsize = 20, figure_padding = 30), axis = (;), plot = (;))

Produces a plot of volume-weighted crystal size distributions. Data are plotted as points while model outputs are shown as lines.

dataset contains the data to be plotted. If providing system, the plot will also include the model output as comparison to the data. ids can be used to specify which datasets are included in the plot.

For mode = :combined all selected datasets will be shown in a single subfigure whereas mode = :separate plots each CSD in its own subfigure.

Figure, Axis and Plot properties can be freely adjusted by providing the keywords figure, axis, and plot with a named tuple, e.g., figure = (fontsize = 20, figure_padding = 30). Informatoin on valid attributes can be found in the Makie.jl Documentation.

plot_attainable_region_dynamic(system, conditions, d43, P; figure = (fontsize = 24, 
    size = (1500, 900), figure_padding = 30), axis = (;), plot = (;))

Produces a plot that compares shows quality of fit by comparing simulated CSDs with the respective dataset on a logarithmic scale. Data are plotted as points while model outputs are shown as lines.

dataset contains the data to be plotted. system contains the kinetic para, the plot will also include the model output as comparison to the data. ids can be used to specify which datasets are included in the plot.

For mode = :combined all selected datasets will be shown in a single subfigure whereas mode = :separate plots each CSD in its own subfigure.

Figure, Axis and Plot properties can be freely adjusted by providing the keywords figure, axis, and plot with a named tuple, e.g., figure = (fontsize = 20, figure_padding = 30). Information on valid attributes can be found in the Makie.jl Documentation.

plot_fit_quality(system, dataset; mode = :combined, ids = 1:length(dataset.T), 
        figure = (fontsize = 20, figure_padding = 30), 
        axis = (xticks = 0:400:2000, yticks = 0:5:35, aspect = 1, limits = (0, 2000, 15, 35)), 
        plot = (;), 
        axislegend = (position = :rt, framevisible = false))

Produces a plot that compares shows quality of fit by comparing simulated CSDs with the respective dataset on a logarithmic scale. Data are plotted as points while model outputs are shown as lines.

dataset contains the data to be plotted. system contains the kinetic para, the plot will also include the model output as comparison to the data. ids can be used to specify which datasets are included in the plot.

For mode = :combined all selected datasets will be shown in a single subfigure whereas mode = :separate plots each CSD in its own subfigure.

Figure, Axis and Plot properties can be freely adjusted by providing the keywords figure, axis, and plot with a named tuple, e.g., figure = (fontsize = 20, figure_padding = 30). Information on valid attributes can be found in the Makie.jl Documentation.

n, Css, d43 = simulate_CSD(L, τ, T, Cf, system::SystemSpecification)

Calculates the CSD n (discretized at sizes L), the steady state concentration Css, and the volum-weighted mean particle size d43 occuring in an MSMPR crystallizer at steady state when it is operated at temperature T, residence time τ, and feed concentration Cf. Kinetics and solubility are specified in system.