From 7107b91a4ec60c31cd34dbf3aa5234b2319476c8 Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Mon, 14 Sep 2026 11:08:12 -0400
Subject: [PATCH 01/18] doc: wip DAE example in manual
---
docs/make.jl | 3 +-
docs/src/manual/mtk.md | 4 +-
.../manual/{nonlinmpc.md => nonlinmpc1.md} | 6 +-
docs/src/manual/nonlinmpc2.md | 132 ++++++++++++++++++
4 files changed, 139 insertions(+), 6 deletions(-)
rename docs/src/manual/{nonlinmpc.md => nonlinmpc1.md} (99%)
create mode 100644 docs/src/manual/nonlinmpc2.md
diff --git a/docs/make.jl b/docs/make.jl
index 4f107581e..0a668960b 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -40,7 +40,8 @@ makedocs(
"Installation" => "manual/installation.md",
"Examples" => [
"Linear Design" => "manual/linmpc.md",
- "Nonlinear Design" => "manual/nonlinmpc.md",
+ "Nonlinear Design (ODE)" => "manual/nonlinmpc1.md",
+ "Nonlinear Design (DAE)" => "manual/nonlinmpc2.md",
"ModelingToolkit" => "manual/mtk.md",
],
],
diff --git a/docs/src/manual/mtk.md b/docs/src/manual/mtk.md
index 6c3b7d1fc..1f942cf2a 100644
--- a/docs/src/manual/mtk.md
+++ b/docs/src/manual/mtk.md
@@ -11,8 +11,8 @@ old_logger = global_logger(); global_logger(errlogger);
## Pendulum Model
-This example integrates the simple pendulum model of the [last section](@ref man_nonlin) in
-[`ModelingToolkit`](https://docs.sciml.ai/ModelingToolkit/stable/) (MTK) framework and
+This example integrates the simple pendulum model of the [penultimate section](@ref man_nonlin)
+in [`ModelingToolkit`](https://docs.sciml.ai/ModelingToolkit/stable/) (MTK) framework and
extracts appropriate `f!` and `h!` functions to construct a [`NonLinModel`](@ref). An
[`NonLinMPC`](@ref) is designed from this model and simulated to reproduce the results of
the last section.
diff --git a/docs/src/manual/nonlinmpc.md b/docs/src/manual/nonlinmpc1.md
similarity index 99%
rename from docs/src/manual/nonlinmpc.md
rename to docs/src/manual/nonlinmpc1.md
index 5ed7e2026..9c6aad088 100644
--- a/docs/src/manual/nonlinmpc.md
+++ b/docs/src/manual/nonlinmpc1.md
@@ -1,7 +1,7 @@
-# [Manual: Nonlinear Design](@id man_nonlin)
+# [Manual: Nonlinear Design (ODE)](@id man_nonlin)
```@contents
-Pages = ["nonlinmpc.md"]
+Pages = ["nonlinmpc1.md"]
```
## Nonlinear Model
@@ -38,7 +38,7 @@ in which ``g`` is the gravitational acceleration in m/s², ``L``, the pendulum l
the end of the pendulum in kg, all bundled in the parameter vector ``\mathbf{p} =
[\begin{smallmatrix} g & L & K & m \end{smallmatrix}]'``. The [`NonLinModel`](@ref)
constructor assumes by default that the state function `f` is continuous in time, that is,
-an ordinary differential equation system (like here):
+an ordinary differential equation (ODE) system (like here):
```@codeblocks
line_counter = :continue
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
new file mode 100644
index 000000000..62b64efc4
--- /dev/null
+++ b/docs/src/manual/nonlinmpc2.md
@@ -0,0 +1,132 @@
+# [Manual: Nonlinear Design (DAE)](@id man_ade)
+
+```@contents
+Pages = ["nonlinmpc2.md"]
+```
+
+The solution must remain electrically neutral, i.e. the sum of the charges of all ions must
+equal zero. Here, the ions are:
+
+- Hydrogen ``[\mathrm{H}^+]``
+- Sodium ``[\mathrm{Na}^+]``
+- Hydroxide ``[\mathrm{OH}^-]``
+- Acetate ``[\mathrm{Ac}^-]``
+
+The charge balance leads to:
+
+```math
+0 = [\mathrm{H}^+] + [\mathrm{Na}^+] - [\mathrm{OH}^-] - [\mathrm{Ac}^-]
+```
+
+The water dissociation constant and the weak acid equilibrium constant are respectively
+defined by:
+
+```math
+\begin{aligned}
+K_w &= [\mathrm{H}^+][\mathrm{OH}^-] \\
+K_a &= \frac{[\mathrm{H}^+][\mathrm{Ac}^-]}{[\mathrm{H}\mathrm{Ac}]}
+\end{aligned}
+```
+
+We have one algebraic variable and two states, respectively denoted with:
+
+```math
+\begin{aligned}
+a_H &= [\mathrm{H}^+] \\
+c_A &= [\mathrm{H}\mathrm{Ac}] + [\mathrm{Ac}^-] \\
+c_B &= [\mathrm{Na}^+]
+\end{aligned}
+```
+
+in which ``c_A`` represents the total concentration of the acid species in the reactor
+(mol/L), composed of an undissociated acid ``\mathrm{H}\mathrm{Ac}`` and the acetate ion
+``\mathrm{Ac}^-``. Substituting the constants, algebraic and state variables into the charge
+balance leads the algebraic equation:
+
+```math
+0 = a_H + c_B - \frac{K_w}{a_H} - \frac{K_a c_A}{K_a + a_H}
+```
+
+The code is:
+
+```@example 1
+using ModelPredictiveControl
+
+# Process parameters
+V = 1000.0 # Reactor volume [L]
+c_Ain = 0.1 # Feed concentration of weak acid [mol/L]
+c_Bin = 0.1 # Feed concentration of strong base [mol/L]
+Kw = 1.0e-14 # Water dissociation constant [mol^2/L^2]
+Ka = 1.75e-5 # Acid dissociation constant [mol/L]
+
+function fq!(ẋ, res, x, a, u, d, p)
+ c_Ain, c_Bin, Kw, Ka, V = p
+ q_Ain, q_Bin = d[1], u[1] # [L/min], [L/min]
+ c_A, c_B = x[1], x[2] # [mol/L], [mol/L]
+ a_H = a[1] # [mol/L]
+ q_out = q_Ain + q_Bin # overflow weir
+ c_Aout = c_A # perfect mixing assumption
+ c_Bout = c_B # perfect mixing assumption
+ ẋ[1] = (q_Ain * c_Ain - q_out * c_Aout) / V
+ ẋ[2] = (q_Bin * c_Bin - q_out * c_Bout) / V
+ res[1] = a_H + c_B - (Kw / a_H) - (Ka * c_A / (Ka + a_H))
+ return nothing
+end
+
+function h!(y, _, a, _ , _ )
+ a_H = a[1]
+ y[1] = try
+ -log10(a_H) # y = pH
+ catch myerror
+ if myerror isa DomainError
+ NaN
+ else
+ rethrow()
+ end
+ end
+ return nothing
+end
+
+Ts = 15.0 # Sample time [min]
+nu = 1 # Manipulated input: base flow rate q_B
+nx = 2 # Differential states: total acid W_a, total base W_b
+na = 1 # Algebraic variable: Hydrogen ion concentration [H+]
+ny = 1 # Measured output: pH
+nd = 1 # Measured disturbance: acid flow rate q_A
+p = [c_Ain, c_Bin, Kw, Ka, V]
+
+model = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny, nd; p)
+model = setname!(model, u=["\$q_B\$"], x=["\$c_a\$", "\$c_b\$"], y=["\$\\mathrm{pH}\$"], d=["\$q_A\$"])
+
+model.a0 .= 1e-5 # initial guess for algebraic variable (positive value because of log10)
+model.Z .= 1e-5 # initial guess for collocation points (positive value because of log10)
+
+u = [10.0]
+d = [10.0]
+x_0 = [0.05556, 0.04444]
+#res = sim!(model, 50, u, d; x_0)
+N = 100
+Y_data, U_data, D_data, X_data = zeros(ny, N), zeros(nu, N), zeros(nd, N), zeros(nx, N)
+x = x_0
+let x=x, u=u, d=d
+ setstate!(model, x)
+ for i=1:N
+ #@show x
+ d = [10.0]
+ y = model(d)
+ u = i < N/2 ? [10.0] : [9.5]
+ Y_data[:, i] = y
+ U_data[:, i] = u
+ D_data[:, i] = d
+ X_data[:, i] = x
+ x = updatestate!(model, u, d)
+ end
+end
+res = SimResult(model, U_data, Y_data, D_data; X_data)
+
+using Plots
+theme(:default)
+#theme(:dark)
+default(fontfamily="Computer Modern"); scalefontsizes(1.1)
+plot(res, plotx=true, plotd=false)
+```
From 90057457fe4b13c69e3608d2108692fdd3e01b07 Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Mon, 14 Sep 2026 11:16:42 -0400
Subject: [PATCH 02/18] doc: minor correction
---
docs/src/manual/nonlinmpc2.md | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index 62b64efc4..34642e5d5 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -28,7 +28,7 @@ K_a &= \frac{[\mathrm{H}^+][\mathrm{Ac}^-]}{[\mathrm{H}\mathrm{Ac}]}
\end{aligned}
```
-We have one algebraic variable and two states, respectively denoted with:
+We respectively denote the algebraic variable and the two states with:
```math
\begin{aligned}
From 458e7a69fac52d3f20587595658e1d5c6d51b413 Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 11:33:24 -0400
Subject: [PATCH 03/18] doc: details on the pH model
---
docs/src/index.md | 3 +-
docs/src/manual/nonlinmpc1.md | 4 +-
docs/src/manual/nonlinmpc2.md | 89 +++++++++++++++++++++++++++--------
3 files changed, 74 insertions(+), 22 deletions(-)
diff --git a/docs/src/index.md b/docs/src/index.md
index 7b7c5e820..a45da5459 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -33,7 +33,8 @@ Depth = 2
Pages = [
joinpath("manual", "installation.md"),
joinpath("manual", "linmpc.md"),
- joinpath("manual", "nonlinmpc.md"),
+ joinpath("manual", "nonlinmpc1.md"),
+ joinpath("manual", "nonlinmpc2.md"),
joinpath("manual", "mtk.md")
]
```
diff --git a/docs/src/manual/nonlinmpc1.md b/docs/src/manual/nonlinmpc1.md
index 9c6aad088..bdbec4107 100644
--- a/docs/src/manual/nonlinmpc1.md
+++ b/docs/src/manual/nonlinmpc1.md
@@ -24,7 +24,7 @@ The following figure presents the system:
border:20px solid white; display: block; margin-left: auto; margin-right: auto;"/>
```
-The plant model is nonlinear:
+The plant model is nonlinear ordinary differential equation (ODE) system:
```math
\begin{aligned}
@@ -38,7 +38,7 @@ in which ``g`` is the gravitational acceleration in m/s², ``L``, the pendulum l
the end of the pendulum in kg, all bundled in the parameter vector ``\mathbf{p} =
[\begin{smallmatrix} g & L & K & m \end{smallmatrix}]'``. The [`NonLinModel`](@ref)
constructor assumes by default that the state function `f` is continuous in time, that is,
-an ordinary differential equation (ODE) system (like here):
+an ODE system (like here):
```@codeblocks
line_counter = :continue
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index 34642e5d5..9bbe874d4 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -4,6 +4,23 @@
Pages = ["nonlinmpc2.md"]
```
+## Nonlinear Model
+
+In this example, the goal is to control the pH of a solution in a continuously stirred
+tank reactor (CSTR) for neutralization. The manipulated input is the inlet flow rate of
+a strong base in L/min, while the inlet flow rate of a weak acid, also in L/min, is a
+measured disturbance:
+
+```math
+\begin{aligned}
+ \mathbf{u} &= q_{Bin} \\
+ \mathbf{y} &= \mathrm{pH} \\
+ \mathbf{d} &= q_{Ain}
+\end{aligned}
+```
+
+### Instantaneous Charge Balance
+
The solution must remain electrically neutral, i.e. the sum of the charges of all ions must
equal zero. Here, the ions are:
@@ -12,7 +29,7 @@ equal zero. Here, the ions are:
- Hydroxide ``[\mathrm{OH}^-]``
- Acetate ``[\mathrm{Ac}^-]``
-The charge balance leads to:
+The steady-state charge balance leads to:
```math
0 = [\mathrm{H}^+] + [\mathrm{Na}^+] - [\mathrm{OH}^-] - [\mathrm{Ac}^-]
@@ -32,9 +49,9 @@ We respectively denote the algebraic variable and the two states with:
```math
\begin{aligned}
-a_H &= [\mathrm{H}^+] \\
-c_A &= [\mathrm{H}\mathrm{Ac}] + [\mathrm{Ac}^-] \\
-c_B &= [\mathrm{Na}^+]
+ a_H &= [\mathrm{H}^+] \\
+ c_A &= [\mathrm{H}\mathrm{Ac}] + [\mathrm{Ac}^-] \\
+ c_B &= [\mathrm{Na}^+]
\end{aligned}
```
@@ -47,6 +64,46 @@ balance leads the algebraic equation:
0 = a_H + c_B - \frac{K_w}{a_H} - \frac{K_a c_A}{K_a + a_H}
```
+The pH is computed with:
+
+```math
+\mathrm{pH} = -10 \log_{10}(a_H)
+```
+
+!!! details "Reduction to an ODE"
+ The algebraic equation can be further manipulated to produce this cubic expression:
+ ```math
+ 0 = a_H^3 + (c_B + K_a) a_H^2 + \big(K_a(c_A + c_B) + K_w \big) a_H - K_w K_a
+ ```
+ We could extract the positive root of this expression inside the output function `h!`
+ to transform the system to an ODE, effectively avoiding the increased complexity of
+ DAEs. The tutorial will still treat the system as a DAE to illustrate its API.
+
+### Mass Balance
+
+Applying a mass balance on the weak acid ``A`` and the strong base ``B`` invariants leads
+to the differential equations:
+
+```math
+\begin{aligned}
+ \dot{c}_A(t) &= \frac{1}{V}(q_{Ain} c_{Ain} - q_{out} c_{Aout}) \\
+ \dot{c}_B(t) &= \frac{1}{V}(q_{Bin} c_{Bin} - q_{out} c_{Bout})
+\end{aligned}
+```
+
+in which the concentrations ``c`` are in mol/L, the tank volume ``V`` in L and the
+volumetric flow rates ``q`` in L/min. By assuming a perfectly mixed reactor and
+an overflow weir to draw the neutralized solution, the following relations evaluate
+the outlet terms:
+
+```math
+\begin{aligned}
+ c_{Aout} &= c_{A} \\
+ c_{Bout} &= c_{B} \\
+ q_{out} &= q_{Ain} + q_{Bin}
+\end{aligned}
+```
+
The code is:
```@example 1
@@ -64,9 +121,9 @@ function fq!(ẋ, res, x, a, u, d, p)
q_Ain, q_Bin = d[1], u[1] # [L/min], [L/min]
c_A, c_B = x[1], x[2] # [mol/L], [mol/L]
a_H = a[1] # [mol/L]
- q_out = q_Ain + q_Bin # overflow weir
- c_Aout = c_A # perfect mixing assumption
- c_Bout = c_B # perfect mixing assumption
+ q_out = q_Ain + q_Bin # [L/min]
+ c_Aout = c_A # [mol/L]
+ c_Bout = c_B # [mol/L]
ẋ[1] = (q_Ain * c_Ain - q_out * c_Aout) / V
ẋ[2] = (q_Bin * c_Bin - q_out * c_Bout) / V
res[1] = a_H + c_B - (Kw / a_H) - (Ka * c_A / (Ka + a_H))
@@ -88,23 +145,17 @@ function h!(y, _, a, _ , _ )
end
Ts = 15.0 # Sample time [min]
-nu = 1 # Manipulated input: base flow rate q_B
-nx = 2 # Differential states: total acid W_a, total base W_b
-na = 1 # Algebraic variable: Hydrogen ion concentration [H+]
-ny = 1 # Measured output: pH
-nd = 1 # Measured disturbance: acid flow rate q_A
+nu, nx, na, ny, nd = 1, 2, 1, 1, 1
p = [c_Ain, c_Bin, Kw, Ka, V]
-model = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny, nd; p)
-model = setname!(model, u=["\$q_B\$"], x=["\$c_a\$", "\$c_b\$"], y=["\$\\mathrm{pH}\$"], d=["\$q_A\$"])
-
-model.a0 .= 1e-5 # initial guess for algebraic variable (positive value because of log10)
-model.Z .= 1e-5 # initial guess for collocation points (positive value because of log10)
+model = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny, nd; p, as_0=[1e-5])
+vu, vd = ["\$q_B\$ (L/min)"], ["\$q_A\$ (L/min)"]
+vx, vy = ["\$c_a\$ (mol/L)", "\$c_b\$ (mol/L)"], ["\$\\mathrm{pH}\$"]
+model = setname!(model, u=vu, x=vx, y=vy, d=vd)
u = [10.0]
d = [10.0]
-x_0 = [0.05556, 0.04444]
-#res = sim!(model, 50, u, d; x_0)
+x_0 = [0.055, 0.045]
N = 100
Y_data, U_data, D_data, X_data = zeros(ny, N), zeros(nu, N), zeros(nd, N), zeros(nx, N)
x = x_0
From 720094e0590aafef8d4de7f828132d33b2ea387b Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 11:53:24 -0400
Subject: [PATCH 04/18] doc: minor detail
---
docs/src/manual/nonlinmpc2.md | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index 9bbe874d4..a5fb279a1 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -75,9 +75,9 @@ The pH is computed with:
```math
0 = a_H^3 + (c_B + K_a) a_H^2 + \big(K_a(c_A + c_B) + K_w \big) a_H - K_w K_a
```
- We could extract the positive root of this expression inside the output function `h!`
- to transform the system to an ODE, effectively avoiding the increased complexity of
- DAEs. The tutorial will still treat the system as a DAE to illustrate its API.
+ We could extract the positive real root of this expression inside the output function
+ `h!` to transform the system to an ODE, effectively avoiding the increased complexity
+ of DAEs. The tutorial will still treat the system as a DAE to illustrate its API.
### Mass Balance
From 2046853f2775d642a69eb052015b221fa090da2d Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 11:56:43 -0400
Subject: [PATCH 05/18] doc: minor detail
---
docs/src/manual/nonlinmpc1.md | 2 +-
docs/src/manual/nonlinmpc2.md | 2 +-
2 files changed, 2 insertions(+), 2 deletions(-)
diff --git a/docs/src/manual/nonlinmpc1.md b/docs/src/manual/nonlinmpc1.md
index bdbec4107..730164fae 100644
--- a/docs/src/manual/nonlinmpc1.md
+++ b/docs/src/manual/nonlinmpc1.md
@@ -4,7 +4,7 @@
Pages = ["nonlinmpc1.md"]
```
-## Nonlinear Model
+## Nonlinear Model (ODE)
In this example, the goal is to control the angular position ``θ`` of a pendulum
attached to a motor. Knowing that the manipulated input is the motor torque ``τ`` in Nm, the
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index a5fb279a1..8f59b8d94 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -4,7 +4,7 @@
Pages = ["nonlinmpc2.md"]
```
-## Nonlinear Model
+## Nonlinear Model (DAE)
In this example, the goal is to control the pH of a solution in a continuously stirred
tank reactor (CSTR) for neutralization. The manipulated input is the inlet flow rate of
From ff84ffd17b877584fa4e5791c79bcaae7cfd5db6 Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 15:11:47 -0400
Subject: [PATCH 06/18] doc: added pH neutralization figure
---
docs/src/assets/ph_neutralization.svg | 109 ++++++++++++++++++++++++++
docs/src/manual/nonlinmpc2.md | 10 ++-
2 files changed, 118 insertions(+), 1 deletion(-)
create mode 100644 docs/src/assets/ph_neutralization.svg
diff --git a/docs/src/assets/ph_neutralization.svg b/docs/src/assets/ph_neutralization.svg
new file mode 100644
index 000000000..9da1422ac
--- /dev/null
+++ b/docs/src/assets/ph_neutralization.svg
@@ -0,0 +1,109 @@
+
+
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index 8f59b8d94..68f8c1ed0 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -19,10 +19,18 @@ measured disturbance:
\end{aligned}
```
+The following figure depicts the pH neutralization process:
+
+```@raw html
+
+```
+
### Instantaneous Charge Balance
The solution must remain electrically neutral, i.e. the sum of the charges of all ions must
-equal zero. Here, the ions are:
+equal zero. The ions in this case study are:
- Hydrogen ``[\mathrm{H}^+]``
- Sodium ``[\mathrm{Na}^+]``
From 729f223b3cee899fc9eeacda25125fad70ef1d8e Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 15:19:17 -0400
Subject: [PATCH 07/18] doc: simpler `h!` function
---
docs/src/manual/nonlinmpc2.md | 11 ++++-------
1 file changed, 4 insertions(+), 7 deletions(-)
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index 68f8c1ed0..ebad389f2 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -140,15 +140,12 @@ end
function h!(y, _, a, _ , _ )
a_H = a[1]
- y[1] = try
- -log10(a_H) # y = pH
+ pH = try
+ -log10(a_H)
catch myerror
- if myerror isa DomainError
- NaN
- else
- rethrow()
- end
+ myerror isa DomainError ? NaN : rethrow()
end
+ y[1] = pH
return nothing
end
From 4d06ba25efac2ef07664c5f2f3326e8d41f1232b Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 15:39:19 -0400
Subject: [PATCH 08/18] doc: minor change
---
docs/src/manual/nonlinmpc2.md | 11 +++++------
1 file changed, 5 insertions(+), 6 deletions(-)
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index ebad389f2..2f991072f 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -117,12 +117,11 @@ The code is:
```@example 1
using ModelPredictiveControl
-# Process parameters
-V = 1000.0 # Reactor volume [L]
-c_Ain = 0.1 # Feed concentration of weak acid [mol/L]
-c_Bin = 0.1 # Feed concentration of strong base [mol/L]
-Kw = 1.0e-14 # Water dissociation constant [mol^2/L^2]
-Ka = 1.75e-5 # Acid dissociation constant [mol/L]
+V = 1000.0 # reactor volume [L]
+c_Ain = 0.1 # feed concentration of weak acid [mol/L]
+c_Bin = 0.1 # feed concentration of strong base [mol/L]
+Kw = 1.0e-14 # water dissociation constant [mol^2/L^2]
+Ka = 1.75e-5 # acid dissociation constant [mol/L]
function fq!(ẋ, res, x, a, u, d, p)
c_Ain, c_Bin, Kw, Ka, V = p
From 98b093d79cd03c0e8033aba821ee04534342a166 Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 16:10:24 -0400
Subject: [PATCH 09/18] doc: plots in hours for pH neutralization
---
docs/src/manual/nonlinmpc2.md | 21 +++++++++++----------
1 file changed, 11 insertions(+), 10 deletions(-)
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index 2f991072f..8b3d27f0f 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -94,15 +94,15 @@ to the differential equations:
```math
\begin{aligned}
- \dot{c}_A(t) &= \frac{1}{V}(q_{Ain} c_{Ain} - q_{out} c_{Aout}) \\
- \dot{c}_B(t) &= \frac{1}{V}(q_{Bin} c_{Bin} - q_{out} c_{Bout})
+ \dot{c}_A(t) &= \frac{60}{V}(q_{Ain} c_{Ain} - q_{out} c_{Aout}) \\
+ \dot{c}_B(t) &= \frac{60}{V}(q_{Bin} c_{Bin} - q_{out} c_{Bout})
\end{aligned}
```
in which the concentrations ``c`` are in mol/L, the tank volume ``V`` in L and the
-volumetric flow rates ``q`` in L/min. By assuming a perfectly mixed reactor and
-an overflow weir to draw the neutralized solution, the following relations evaluate
-the outlet terms:
+volumetric flow rates ``q`` in L/min. The accumulation terms ``\dot{c}`` are in mol/(L h)
+because of the ``60`` factor. By assuming a perfectly mixed reactor and an overflow weir to
+draw the neutralized solution, the following relations evaluate the outlet terms:
```math
\begin{aligned}
@@ -131,8 +131,8 @@ function fq!(ẋ, res, x, a, u, d, p)
q_out = q_Ain + q_Bin # [L/min]
c_Aout = c_A # [mol/L]
c_Bout = c_B # [mol/L]
- ẋ[1] = (q_Ain * c_Ain - q_out * c_Aout) / V
- ẋ[2] = (q_Bin * c_Bin - q_out * c_Bout) / V
+ ẋ[1] = (60/V)*(q_Ain * c_Ain - q_out * c_Aout)
+ ẋ[2] = (60/V)*(q_Bin * c_Bin - q_out * c_Bout)
res[1] = a_H + c_B - (Kw / a_H) - (Ka * c_A / (Ka + a_H))
return nothing
end
@@ -148,7 +148,7 @@ function h!(y, _, a, _ , _ )
return nothing
end
-Ts = 15.0 # Sample time [min]
+Ts = 0.5 # Sample time [h]
nu, nx, na, ny, nd = 1, 2, 1, 1, 1
p = [c_Ain, c_Bin, Kw, Ka, V]
@@ -160,7 +160,7 @@ model = setname!(model, u=vu, x=vx, y=vy, d=vd)
u = [10.0]
d = [10.0]
x_0 = [0.055, 0.045]
-N = 100
+N = 50
Y_data, U_data, D_data, X_data = zeros(ny, N), zeros(nu, N), zeros(nd, N), zeros(nx, N)
x = x_0
let x=x, u=u, d=d
@@ -183,5 +183,6 @@ using Plots
theme(:default)
#theme(:dark)
default(fontfamily="Computer Modern"); scalefontsizes(1.1)
-plot(res, plotx=true, plotd=false)
+p = plot(res, plotx=true, plotd=false, xlabel="Time (h)")
+xlabel!(p[3], "")
```
From c55e8b39323485ddbff015ce095ad1375a4026a8 Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 16:11:57 -0400
Subject: [PATCH 10/18] doc: correct names
---
docs/src/manual/nonlinmpc2.md | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index 8b3d27f0f..b01b4260a 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -154,7 +154,7 @@ p = [c_Ain, c_Bin, Kw, Ka, V]
model = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny, nd; p, as_0=[1e-5])
vu, vd = ["\$q_B\$ (L/min)"], ["\$q_A\$ (L/min)"]
-vx, vy = ["\$c_a\$ (mol/L)", "\$c_b\$ (mol/L)"], ["\$\\mathrm{pH}\$"]
+vx, vy = ["\$c_A\$ (mol/L)", "\$c_B\$ (mol/L)"], ["\$\\mathrm{pH}\$"]
model = setname!(model, u=vu, x=vx, y=vy, d=vd)
u = [10.0]
From 9a48be98eb0ab1bb265ba671c0af7c7a9b2703bb Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 16:28:53 -0400
Subject: [PATCH 11/18] doc: new litmus colors to be more intuitive
---
docs/src/assets/ph_neutralization.svg | 34 +++++++++++++--------------
1 file changed, 17 insertions(+), 17 deletions(-)
diff --git a/docs/src/assets/ph_neutralization.svg b/docs/src/assets/ph_neutralization.svg
index 9da1422ac..ac0246129 100644
--- a/docs/src/assets/ph_neutralization.svg
+++ b/docs/src/assets/ph_neutralization.svg
@@ -46,10 +46,10 @@
-
-
-
-
+
+
+
+
@@ -68,34 +68,34 @@
-
-
-
-
-
+
+
+
+
+
-
+
-
+
-
+
-
-
+
+
-
+
-
-
+
+
From fb896772ce39ba5038806513f41cc556c9b6c69e Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 16:29:52 -0400
Subject: [PATCH 12/18] doc: minor change
---
docs/src/manual/nonlinmpc2.md | 1 -
1 file changed, 1 deletion(-)
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index b01b4260a..55950c259 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -166,7 +166,6 @@ x = x_0
let x=x, u=u, d=d
setstate!(model, x)
for i=1:N
- #@show x
d = [10.0]
y = model(d)
u = i < N/2 ? [10.0] : [9.5]
From 607f3b5998df79388a4121a473656118f2553dbf Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 16:44:44 -0400
Subject: [PATCH 13/18] doc: minor tweaks
---
docs/src/manual/nonlinmpc2.md | 17 +++++++++--------
1 file changed, 9 insertions(+), 8 deletions(-)
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index 55950c259..1c15535b7 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -159,16 +159,16 @@ model = setname!(model, u=vu, x=vx, y=vy, d=vd)
u = [10.0]
d = [10.0]
-x_0 = [0.055, 0.045]
-N = 50
+x_0 = [0.051, 0.049]
+N = 61
Y_data, U_data, D_data, X_data = zeros(ny, N), zeros(nu, N), zeros(nd, N), zeros(nx, N)
x = x_0
let x=x, u=u, d=d
setstate!(model, x)
for i=1:N
- d = [10.0]
+ d = i ≤ 2N÷3 ? [10.0] : [9.8]
y = model(d)
- u = i < N/2 ? [10.0] : [9.5]
+ u = i ≤ N÷3 ? [10.0] : [9.7]
Y_data[:, i] = y
U_data[:, i] = u
D_data[:, i] = d
@@ -179,9 +179,10 @@ end
res = SimResult(model, U_data, Y_data, D_data; X_data)
using Plots
-theme(:default)
-#theme(:dark)
+#theme(:default)
+theme(:dark)
default(fontfamily="Computer Modern"); scalefontsizes(1.1)
-p = plot(res, plotx=true, plotd=false, xlabel="Time (h)")
-xlabel!(p[3], "")
+#p = plot(res, plotx=true, plotd=false, xlabel="Time (h)")
+#xlabel!(p[3], "")
+p = plot(res, plotd=true, xlabel="Time (h)")
```
From 6bba28133d3cd4451819eafcd07a1784e3f5b9de Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 16:46:26 -0400
Subject: [PATCH 14/18] doc: minor tweaks
---
docs/src/manual/nonlinmpc2.md | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index 1c15535b7..fdbba6c59 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -101,7 +101,7 @@ to the differential equations:
in which the concentrations ``c`` are in mol/L, the tank volume ``V`` in L and the
volumetric flow rates ``q`` in L/min. The accumulation terms ``\dot{c}`` are in mol/(L h)
-because of the ``60`` factor. By assuming a perfectly mixed reactor and an overflow weir to
+because of the 60 factors. By assuming a perfectly mixed reactor and an overflow weir to
draw the neutralized solution, the following relations evaluate the outlet terms:
```math
From 54a35415052d3801209ca3dbc83d492e5f13b927 Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 16:56:42 -0400
Subject: [PATCH 15/18] doc: minor correction
---
docs/src/manual/nonlinmpc2.md | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index fdbba6c59..be0d05f5e 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -37,7 +37,7 @@ equal zero. The ions in this case study are:
- Hydroxide ``[\mathrm{OH}^-]``
- Acetate ``[\mathrm{Ac}^-]``
-The steady-state charge balance leads to:
+The instantaneous charge balance leads to:
```math
0 = [\mathrm{H}^+] + [\mathrm{Na}^+] - [\mathrm{OH}^-] - [\mathrm{Ac}^-]
From 630168e7e72ac01409c794dd73432f06c3637ce6 Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Tue, 15 Sep 2026 17:12:33 -0400
Subject: [PATCH 16/18] doc: minor corrections
---
docs/src/manual/nonlinmpc2.md | 9 +++++----
1 file changed, 5 insertions(+), 4 deletions(-)
diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index be0d05f5e..068ec43f3 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -1,4 +1,4 @@
-# [Manual: Nonlinear Design (DAE)](@id man_ade)
+# [Manual: Nonlinear Design (DAE)](@id man_dae)
```@contents
Pages = ["nonlinmpc2.md"]
@@ -19,7 +19,8 @@ measured disturbance:
\end{aligned}
```
-The following figure depicts the pH neutralization process:
+An overflow weir draws the neutralized solution, effectively keeping a constant volume of
+solution inside the tank. The following figure depicts the pH neutralization process:
```@raw html
Date: Wed, 16 Sep 2026 18:27:09 -0400
Subject: [PATCH 17/18] doc: green solution in tank (neutral Litmus color)
---
docs/src/assets/ph_neutralization.svg | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/docs/src/assets/ph_neutralization.svg b/docs/src/assets/ph_neutralization.svg
index ac0246129..29c50f9b3 100644
--- a/docs/src/assets/ph_neutralization.svg
+++ b/docs/src/assets/ph_neutralization.svg
@@ -46,10 +46,10 @@
-
-
-
-
+
+
+
+
From b2a189a0dcb70a42b4dac1e5997be9feafd0c65f Mon Sep 17 00:00:00 2001
From: franckgaga
Date: Thu, 17 Sep 2026 08:33:25 -0400
Subject: [PATCH 18/18] doc: permutate blue and red line on pH tank
---
docs/src/assets/ph_neutralization.svg | 90 +++++++++++++--------------
1 file changed, 45 insertions(+), 45 deletions(-)
diff --git a/docs/src/assets/ph_neutralization.svg b/docs/src/assets/ph_neutralization.svg
index 29c50f9b3..08d106345 100644
--- a/docs/src/assets/ph_neutralization.svg
+++ b/docs/src/assets/ph_neutralization.svg
@@ -37,73 +37,73 @@