Chemical equilibrium

Author: Azzaare

src/ConstraintModels.jl

@@ -19,6 +19,7 @@ export sudoku
 export chemical_equilibrium
 
 include("assignment.jl")
+include("chemical_equilibrium.jl")
 include("cut.jl")
 include("golomb.jl")
 include("magic_square.jl")

src/chemical_equilibrium.jl

@@ -8,10 +8,15 @@ function chemical_equilibrium(A, B, C)
     @variable(model, 0 ≤ X[1:n] ≤ maximum(B))
 
     # mass_conservation function
-    conserve = i -> (x -> sum(A[:, i] .* x) == B[i])
+    conserve = i -> (x ->
+        begin
+            δ = abs(sum(A[:, i] .* x) - B[i])
+            return δ ≤ 1.e-6 ? 0. : δ
+        end
+    )
 
     for i in 1:m
-        @constraint(model, X in Predicate(conserve(i)))
+        @constraint(model, X in Error(conserve(i)))
     end
 
     # computes the total energy freed by the reaction

test/JuMP.jl

@@ -69,3 +69,11 @@ end
 
     @info "JuMP: basic opt" value(x) value(y) (12*value(x)+20*value(y))
 end
+
+@testset "JuMP: Chemical equilibrium" begin
+    m, X = chemical_equilibrium(atoms_compounds, elements_weights, standard_free_energy)
+    # set_optimizer_attribute(m, "iteration", 10000)
+    # set_time_limit_sec(m, 120.0)
+    optimize!(m)
+    @info "JuMP: $compounds_names ⟺ $mixture_name" value.(X)
+end

test/instances.jl

@@ -69,4 +69,8 @@ standard_free_energy = [
     -22.179
 ]
 
-compounds_names = ["H", "H₂", "H₂O", "N", "N₂", "NH", "NO", "O", "O₂", "OH"]
+compounds_names = "x₁⋅H + x₂⋅H₂ + x₃⋅H₂O + x₄⋅N + x₅⋅N₂ + x₆⋅NH + x₇⋅NO + x₈⋅O + x₉⋅O₂ + x₁₀⋅OH"
+
+mixture_name =  "½⋅N₂H₄ + ½⋅O₂"
+
+equation = compounds_names * " = " * mixture_name