How can I solve this 4-equation NonlinearSystem? It gives `retcode: Unstable`

Specific Domains Modelling & Simulations
,
1

I am trying to solve a system of equations with NonlinearSolve.jl and ModelingToolkit.jl.

Specifically a simplified version of Stewart’s model for Acid-Base disturbances in blood (Independence of Variables in Stewart's Model of the Acid-Base Chemistry of the Blood Plasma - ScienceDirect).

The following code returns a meaningless result (a negative H+ concentration), and retcode: Unstable.

using NonlinearSolve
using ModelingToolkit

# constants
Kc = 2.46e-11 
Ka = 3e-7

@variables HA A H HCO3
@parameters pCO2 Atot SID  

# balance equations
balance_eqs = [
	Atot ~ HA + A
	SID ~ HCO3 + A
]

# mass action equations
action_eqs = [
	H ~ Ka * (HA/A)
	H ~ Kc * (pCO2/HCO3)	
]

stewart_eqs = [balance_eqs; action_eqs]

@mtkbuild ns_stewart = NonlinearSystem(stewart_eqs)

# Initial values (guesses)
guesses_stewart = (
	HCO3 => 24.9,
	H => 3.95e-8,
	A => 12.1,
	HA => 5.4
)

ps = (
    Atot => 20,
    SID => 20,
    pCO2 => 5
)

prob = NonlinearProblem(ns_stewart, guesses_stewart, ps)

sol = solve(prob)
2 
julia> full_equations(ns_stewart)
1-element Vector{Equation}:
 0 ~ (-2.46e-11pCO2) / HCO3 + (3.0e-7(Atot + HCO3 - SID)) / (-HCO3 + SID)

julia> unknowns(ns_stewart)
1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 HCO3

From this you can see that the only solution is HCO3 -> infty

2 Likes
3

Thanks! I will try to figure out where the mistake is, and hopefully post a working model.

4

MTK works correct.
Please be more carefully with units of the parameters and initial values you are using. e.g.

  • parameter Kc is using mm-mercury => 1mmHg = 133.322Pa
  • pressure pCO2 = 5.7kPa - you used 5Pa as initial value
  • concentrations should also be checked
1 Like
5

a revised version below:

using NonlinearSolve
using ModelingToolkit

# !!! used UNITS !!!
#  eq/liter
#  Pa

# constants
mmHG = 133.322 # Pa
Kc = 2.46e-11/mmHG
Ka = 3e-7
@variables HA A H HCO3
@parameters pCO2 Atot SID

eqs = [
	H * A ~ Ka * HA
	H * HCO3 ~ Kc * pCO2
	Atot ~ HA + A
	SID ~ HCO3 + A
]

@mtkbuild ns_stewart = NonlinearSystem(eqs)
display(ns_stewart)
println()

# Initial values (guesses)
guesses_stewart = (
	HCO3 => 24.9e-3,
	H => 3.95e-8,
	A => 17.1e-3,
	HA => 2e-3
)

ps = (
    Atot => 19e-3,
    SID => 41.7e-3,
    pCO2 => 5.3e3
)

prob = NonlinearProblem(ns_stewart, guesses_stewart, ps)
sol = solve(prob)

@show sol
@show H = -log10(sol[1])
println("HCO3 = $(sol[2]*1000) meq/liter")

NEW Results:

julia> include("stewart1.jl")
Model ns_stewart with 2 equations
Unknowns (2):
  H
  HCO3
Parameters (3):
  pCO2
  Atot
  SID

sol = [3.9275240698625016e-8, 0.02489948137605699]
H = -(log10(sol[1])) = 7.405881144702646
HCO3 = 24.89948137605699 meq/liter
1 Like
6

Thanks! Reordering the equations fixes it. Also, good catch with the units.

However, I still do not understand why this specification works:

eqs1 = [
	H * A ~ Ka * HA
	H * HCO3 ~ Kc * pCO2
	Atot ~ HA + A
	SID ~ HCO3 + A
]

@mtkbuild mod1 = NonlinearSystem(eqs1)
full_equations(mod1)
#  2-element Vector{Equation}:
#    0 ~ 1.8451568383312584e-10pCO2 - H*HCO3
#    0 ~ -Atot - HCO3 + SID + 3.3333333333333335e6H*(-HCO3 + SID)

while this does not work:

eqs2 = [
	H ~ Ka * HA / A
	H ~ Kc * pCO2 / HCO3
	Atot ~ HA + A
	SID ~ HCO3 + A
]

@mtkbuild mod2 = NonlinearSystem(eqs2)
full_equations(mod2)
# 1-element Vector{Equation}:
#     0 ~ (1.8451568383312584e-10pCO2) / HCO3 + (-3.0e-7(Atot + HCO3 - SID)) / (-HCO3 + SID)

To me, they seem equivalent, but is there something special about how MTK handles / ?

7

My interpretation:
In equ1 I used the equations I found in the article you referred. The variable H here is a part of the LHS expression.
In equ2 the variable H is explicit mentioned in 2 equations. For that reason H was eliminated by MTK.

1 Like
8

In one case H is nonlinearly related to another variable. It’s just a case the tearing algorithm does not handle yet. In fact, it’s hard to reduce because it’s not generally valid to divide like that because in theory that unknown could be zero.

1 Like
9

Thanks for the explanation. I hope I’ll be able to recognize it if I run into a similar problem in the future.