Hi Everbody!
I’m still quite new to Julia, so I hope I am not missing something obvious here…
My task is to solve numerically a system of equations, my working example looks like this:
using Optimization, Zygote, Symbolics
using ForwardDiff
function getChemicalPotential(id::Int, T)
n = ["C(gr)", "H2O", "CH4", "N2", "C2H6", "NH3", "H2", "CO", "CO2", "CH2O2", "C2H4", "C(d)", "HCN"]
c_vector = [
[-5214.9, 13.037, 0.0051, -8.4343e-7, 5.5694e-11, -74.387],
[-9156, 27.216, 0.0091186, -1.2446e-6, 6.847e-11, 28.25],
[-11271, 25.806, 0.02905, -4.7955e-6, 2.9592e-10, 18.67],
[-8463.6, 26.815, 0.0037938, -6.0948e-7, 3.751e-11, 36.65],
[-17885, 42.55, 0.049247, -8.008e-6, 4.7591e-10, -41.251],
[-9535.2, 26.177, 0.019086, -2.9443e-6, 1.6862e-10, 32.5247],
[-7865.8, 26.239, 0.0023634, -1.3255e-7, 1.72229e-12, -20.3209],
[-9086.50436189, 29.0370119629, 0.00245170237304, -2.91117727464e-7, 1.29678792586e-11, 29.8419],
[-13381, 43.301, 0.005694, -6.7725e-7, 3.0335e-11, -38.8709],
[-17807, 48.381, 0.024825, -4.561e-6, 3.0382e-10, -40.997],
[-14897, 37.224, 0.035119, -5.9621e-6, 3.7559e-10, -13.0221],
[-4832, 10.818, 0.0067481, -1.1898e-6, 7.6844e-11, -66.141],
[-10962, 33.517, 0.010353, -1.6675e-6, 1.0168e-10, 4.7979],
]
Hf0_vector = [20.1, -241.826, -74.873, 0.0, -84.6576, -45.898, 0.0, -110.529, -393.522, -378.561, 52.467, 28.2, 135.143]
c = c_vector[id]
Hf0 = Hf0_vector[id]
mu = c[1] - c[6]*T - c[2]*T*log(T-1) - c[3]*T^2 - c[4]*T^(3/2) - c[5]*T^(4/3) + Hf0
return mu
end
function minimizationFunction(x::AbstractVector{Tx}, b::AbstractVector{Tb}) where {Tx, Tb}
α = 0.5
β = 0.38
ε = 1.0
θ = 4350
κ = 9.23
R = 8.314
p = x[end-2]
t = x[end-1]
V = x[end]
mu_over_RT = [getChemicalPotential(i, t)/ R / t for i in 1:13]
xtot = sum(x[1:end-3])
funcs = [mu_over_RT[i] + log(x[i]/xtot) + 1 - x[i]/xtot for i in 1:13]
# BKW equation
kovolumes = [0.0, 246.0, 479.0, 392.0, 832.0, 405.0, 110.0, 390.0, 590.0, 663.0, 700.0, 0.0, 510.0]
xBKW = κ*sum([x[i]*kovolumes[i] for i in 1:13])/(xtot*V * (t + θ)^α)
BKW = p*V/R/t - 1 - xBKW * exp(β*xBKW)
append!(funcs, BKW)
func = sum(abs.(funcs))
return func
end
function solve_equilibrium(b::Vector{Float64})
x0 = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1e5, 1000.0, 1.0]
lower_bounds = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
upper_bounds = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1e10, 20000.0, 1000.0]
lower_constr = [0.0, 0.0, 0.0, 0.0]
upper_constr = [0.0, 0.0, 0.0, 0.0]
alphas = [1.0 0.0 0.0 0.0
0.0 2.0 1.0 0.0
1.0 4.0 0.0 0.0
0.0 0.0 0.0 2.0
2.0 6.0 0.0 0.0
0.0 3.0 0.0 1.0
0.0 2.0 0.0 0.0
1.0 0.0 1.0 0.0
1.0 0.0 2.0 0.0
1.0 2.0 2.0 0.0
2.0 4.0 0.0 0.0
1.0 0.0 0.0 0.0
1.0 1.0 0.0 1.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0]
constraint(res, x, b) = (res .= (alphas' * x) - b)
objective = OptimizationFunction(minimizationFunction, Optimization.AutoForwardDiff(), cons=constraint)
opt_problem = OptimizationProblem(objective, x0, b, lb=lower_bounds, ub=upper_bounds, lcons=lower_constr, ucons=upper_constr)
result = solve(opt_problem, Optimization.LBFGS())
end
# C, H, O, N
b = [5.972, 9.317, 1.917, 1.862]
solution = solve_equilibrium(b)
println("Solution vector (x):")
println(solution)
But now, I get the following error:
LoadError: MethodError: no method matching (::Colon)(::Int64, ::Nothing)
The function `Colon()` exists, but no method is defined for this combination of argument types.
Closest candidates are:
(::Colon)(::T, ::Any, ::T) where T<:Real
@ Base range.jl:50
(::Colon)(::A, ::Any, ::C) where {A<:Real, C<:Real}
@ Base range.jl:10
(::Colon)(::T, ::Any, ::T) where T
@ Base range.jl:49
...
I think the problem is my constraint definition, as I don’t get the error if I leave it out, but I can’t figure out why. Can anybody help me?
Thanks in advance!