using OptimizationOptimisers
using Optimization
loss(θ) = norm((predict(θ) - data_train) .* weights)
optf = OptimizationFunction((x, p) -> loss(x), AutoReverseDiff(true))
optprob = OptimizationProblem(optf, θ)
opt = Adam(0.05, (1.0, 0.985)) # η: learning rate
@time results = Optimization.solve(
optprob, opt, callback = callback, maxiters = 180
)
gives error:
ERROR: Function argument passed to autodiff cannot be proven readonly.
If the the function argument cannot contain derivative data, instead call autodiff(Mode, Const(f), ...)
Original code is here
Full code
include("make_OCFEM.jl")
"""
====================================================================
Building OCFEM (orthogonal collocation on finite element method)
for z discretization with cubic hermite polynomials
====================================================================
"""
begin
Ne = 42 # Number of finite elements
N = 2 # Number of interior collocation points
n_components = 1 # Number of chemical species
n_phases = 2 # 2 phases → 1 liquid(c profile) + 1 solid(q profile)
p_order = 4 # Polynomial order + 1
n_variables = n_components * n_phases * (p_order + 2 * Ne - 2)
xₘᵢₙ = 0.0
xₘₐₓ = 1.0 # z domain limits
h = (xₘₐₓ - xₘᵢₙ) / Ne # finite elements' sizes
end
H, A, B = make_OCFEM(Ne, n_phases, n_components)
# H are the polynomial evaluations at collocation points,
# A first derivative, B second derivative
#Building mass matrix
MM = BitMatrix(make_MM(Ne, n_phases, n_components)) # make mass matrix
"""
===================================================================
Build UDE and solve & visualization
===================================================================
"""
# -------------------Defining chromatography parameters---------------------
begin
Qf = 5.0e-2 # Flow rate (dm^3/min)
d = 0.5 # Column Diameter (dm)
L = 2.0 # Bed Length (dm)
a = pi*d^2/4 # Column area (dm^2)
ϵ = 0.5 # Bed porosity
v = Qf/(a*ϵ) # Interstitial velocity (dm/min)
Pe = 21.095632695978704 # Peclet Number
Da_x = v*L/Pe # Axial dispersion (dm^2/min)
c_in = 5.5 # Feed concentration (mg/L)
k_transf = 0.22 # Mass transfer coefficient (1/min)
k_iso = 1.8 # Isotherm affinity parameter (L/mg)
qmax = 55.54 # Isotherm saturation parameter (mg/L)
q_test = qmax * k_iso * c_in^1.0 / (1.0 + k_iso * c_in^1.0)
# Scale parameter for amount adsorbed in solid phase
# Change exponent according to each isotherm: 1.0 for langmuir, 1.5 for sips in this work
end
#----------Define the derivative matrices stencil and node index-----------------
using LinearAlgebra
begin
# Internal node index
l_idx = 2
u_idx = 2Ne + 1
# Boundary node index
lb_idx = 1
ub_idx = 2Ne + 2
∂x = Array(A * inv(H)) # c(x, t) = H*x and ∂x(c(x, t)) = A*x = A*inv(H)*c(x, t)
∂x² = Array(B * inv(H)) # c(x, t) = H*x and ∂x²(c(x, t)) = B*x = B*inv(H)*c(x, t)
end
#-------------------Importing experimental data-------------------
# Vermeulen’s kinetics and langmuir isotherm
using DelimitedFiles
c_exp_data = readdlm("traindata_improved_quad_lang_2min.csv", ',')
# ------------------Initializing Neural networks----------------------
using Random, Lux, ComponentArrays
nn = Chain(
Dense(2, 22, tanh_fast),
Dense(22, 1)
)
ps, st = Lux.setup(Random.default_rng(), nn)
#---------------------building rhs function for DAE solver---------------------
mutable struct Hybrid{}
L::Float64
h::Float64
v::Float64
ϵ::Float64
c_in::Float64
∂x::Matrix{Float64}
∂x²::Matrix{Float64}
Da_x::Float64
k_iso::Float64
qmax::Float64
q_test::Float64
end
using UnPack
function (Hybrid!::Hybrid)(du, u, p, t)
#Aliasing parameters of chromatography model
@unpack L, h, v, ∂x, ∂x², Da_x, ϵ, c_in, k_iso, qmax, q_test = Hybrid!
#---------------------Mass Transfer and equilibrium -----------------
c = (@view u[lb_idx:ub_idx])
q = u[ub_idx+1:2*ub_idx] / q_test # scaling dependent variables
q_eq = qmax * k_iso * c.^1.0 ./ (1.0 .+ k_iso * c.^1.0) / q_test# Change exponent according to each isotherm
x1x2 = [q_eq q]'
∂x_u = ∂x*u
∂x²_u = ∂x²*u
# Neural network output
nn_out = nn(x1x2, p, st)[1] # p of function Hybrid! stands for parameters of Neural network
nn_internal = nn_out[2:end-1]
#------------------------Mass balance---------------------------
# Concentration profile for internal node. du[2:u_idx] for ∂t(c(x, t))
du[l_idx:u_idx] = - (1 - ϵ) / ϵ * nn_internal .- v/(h*L) * ∂x_u[l_idx:u_idx] .+ Da_x / (h*L)^2 * ∂x²_u[l_idx:u_idx]
# Boundary node equations. left multiply with MM should be 0 (du[1] = du[ub_idx] = 0)
du[1] = Da_x / (h*L) * ∂x_u[1] - v * (c[1] - c_in)
du[ub_idx] = ∂x_u[ub_idx] / (h*L)
# Absorption profile (solid phase)
du[ub_idx+1:2*ub_idx] = nn_out # du[ub_idx+1:end] for ∂t(q(x, t))
nothing
end
using OrdinaryDiffEq
rhs = Hybrid(L, h, v, ϵ, c_in, ∂x, ∂x², Da_x, k_iso, qmax, q_test)
chroma_func = ODEFunction(rhs, mass_matrix = MM)
#-----------------------Build and Solve ODE Problem-------------------------
begin
t0, t_end = first(c_exp_data[:, 1]), last(c_exp_data[: , 1])
tspan = (t0, t_end)
datasize = size(c_exp_data, 1)
t_steps = range(t0, t_end; length = datasize)
ps_ca = ComponentArray(ps)
c0 = 1e-3
u0 = c0 * ones(n_variables)
u0[Int(n_variables/2) + 1:end] .= qmax * k_iso * c0^1.0 / (1.0 + k_iso * c0^1.0)
end
chroma_prob = ODEProblem(chroma_func, u0, tspan, ps_ca)
LinearAlgebra.BLAS.set_num_threads(5)
@time sol_init = Array(solve(chroma_prob, FBDF(autodiff = false), saveat = t_steps))
# FBDF(autodiff = false) also works
using Plots
begin
plot(t_steps, sol_init[Int(n_variables/2), :])
scatter!(c_exp_data[:, 1], c_exp_data[:, 2])
end
"""
===================================================================
Training Neural Network
===================================================================
"""
#-----------------------Training Neural Network----------------------
using SciMLSensitivity
function predict(p)
sensealg = QuadratureAdjoint(autojacvec = ReverseDiffVJP(true))
# Interpolating adjoint works well too
prob_ = remake(chroma_prob; p = p)
s_new = Array(solve(prob_, FBDF(autodiff = false), abstol = 5e-7, reltol = 5e-7, saveat = t_steps, sensealg = sensealg))
s_new[Int(n_variables / 2), :] / c_in
end
# Setting up training data
data_train = c_exp_data[:, 2] / c_in
weights = ones(size(data_train, 1))
#--------------------------Loss function---------------------------
function loss(p)
loss = norm((predict(p) - data_train) .* weights)
return loss
end
#-------------------------testing gradients-----------------------
p = copy(ps_ca)
@time loss(p)
@time predict(p)
using Zygote
@time grad_reverse = Zygote.gradient(loss, p)
#--------------------Maximum Likelihood estimation--------------------
using Optimization, DiffEqFlux
adtype = AutoZygote()
optf = OptimizationFunction((x, p) -> loss(x), adtype)
optprob = OptimizationProblem(optf, p)
iter = 1
callback1 = function(p, l)
global iter
println(l)
println(iter)
iter += 1
l < 1.0e-3
end
using OptimizationOptimisers
opt = Adam(0.05, (1.0, 0.985)) # η: learning rate
@time results = Optimization.solve(
optprob, opt, callback = callback1, maxiters = 180
)
make_OCFEM.jl (6.5 KB)
traindata_improved_quad_lang_2min.csv