This notebook builds a stage-wise (tanks-in-series) equilibrium model of a chromatography column and uses it both to generate synthetic data and to recover the model's parameters from that data.
​
The column is discretized into 50 well-mixed equilibrium stages. For two solutes A and B, each governed by its own linear equilibrium constant (KA, KB), a system of first-order ODEs describing stage-to-stage concentration transfer is assembled and solved with NDSolve for a rectangular pulse injection at the column inlet. The resulting exit-concentration profile from the final stage is sampled over time and corrupted with uniform random noise to produce a synthetic "noisy chromatogram", mimicking real detector data.
​
The noisy data is then fed into NonlinearModelFit (using the interior-point method, with bounds around the true KA/KB values) to re-estimate the equilibrium constants from the simulated measurements, and the fitted model is compared against the noisy data and the noise-free simulation to check how well the original parameters are recovered.
Equilibrium constants
In[]:=
KA=0.9;KB=1.3;
Total volume
In[]:=
V
tot
=Lbed*Pi*Dc*Dc/4.;
50 well-mixed stirred tanks or 50 perfect equilibrium stages
In[]:=
NS=50;
Stage volume
In[]:=
V
stage
=
V
tot
/NS;
Volumetric flow rate
In[]:=
L=0.4;
Void fraction
In[]:=
ϵ=0.4;
Column diameter
In[]:=
Dc=4;
Column length
In[]:=
Lbed=25;
Governing equations
In[]:=
eq[i_]:=D[
CA
i
[t],{t,1}]L(
CA
i-1
[t]-
CA
i
[t])ϵ
V
stage
+(1-ϵ)
V
stage
KA;eq2[i_]:=D[
CB
i
[t],{t,1}]L(
CB
i-1
[t]-
CB
i
[t])ϵ
V
stage
+(1-ϵ)
V
stage
KB;​​eqn=Table[eq[i],{i,1,50}];
In[]:=
initialCond=Table[
CA
i
[0]0,{i,1,50}];​​sys=Flatten[Join[eqn,initialCond]];​​Unknown=Table[
CA
i
,{i,1,50}];​​eqn2=Table[eq2[i],{i,1,50}];​​initialCond2=Table[
CB
i
[0]0,{i,1,50}];​​sys2=Flatten[Join[eqn2,initialCond2]];​​Unkn2=Table[
CB
i
,{i,1,50}];
Rectangular pulse injection at the column inlet
In[]:=
CA
0
[t_]:=
CA
init
(UnitStep[t]-UnitStep[t-10]);​​
CA
init
=0.23;​​
CB
0
[t_]:=
CB
init
(UnitStep[t]-UnitStep[t-10]);​​
CB
init
=0.23;
In[]:=
sol1=NDSolve[sys,Unknown,{t,0,5000}];sol2=NDSolve[sys2,Unkn2,{t,0,3000}];
Generate synthetic noisy chromatogram using exit concentrations of species A and B
In[]:=
datareal1=Table[{i,First[Evaluate[
CA
50
[t]]/.sol1/.t->i]+0.0001RandomReal[{-1,1}]},{i,0,1600,10}];
In[]:=
datareal2=Table[{i,First[Evaluate[
CB
50
[t]]/.sol2/.t->i]+0.0001RandomReal[{-1,1}]},{i,0,1600,10}];
p1=ListPlot{datareal1,datareal2},

Out[]=
CA
50
CB
50
Use the synthetic noisy chromatogram to fetch KA and KB (the equilibrium constants). To do so, we use ParametricNDSolveValue and NonlinearModelFit
In[]:=
Clear[KA,KB]
In[]:=
SYS=Join[sys,sys2];
In[]:=
Length@SYS
Out[]=
200
In[]:=
UNK=Join[Unknown,Unkn2];
In[]:=
Length@UNK
Out[]=
100
In[]:=
pfunall=ParametricNDSolveValue[SYS,UNK,{t,0,2000},{KA,KB}]​​
Out[]=
ParametricFunction
Expression: {
CA
1
,
CA
2
,
CA
3
,
CA
4
,
CA
5
,
CA
6
,
CA
7
,
CA
8
,
CA
9
,
CA
10
,
CA
11
,
CA
12
,
CA
13
,
CA
14
,
CA
15
,
CA
16
,
CA
17
,
CA
18
,
CA
19
,
CA
20
,
CA
21
,
CA
22
,
CA
23
,
CA
24
,
CA
25
,
CA
26
,
CA
27
,
CA
28
,
CA
29
,
CA
30
,
CA
31
,
CA
32
,
CA
33
,
CA
34
,
CA
35
,
CA
36
,
CA
37
,
CA
38
,
CA
39
,
CA
40
,
CA
41
,
CA
42
,
CA
43
,
CA
44
,
CA
45
,
CA
46
,
CA
47
,
CA
48
,
CA
49
,
CA
50
,
CB
1
,
CB
2
,
CB
3
,
CB
4
,
CB
5
,
CB
6
,
CB
7
,
CB
8
,
CB
9
,
CB
10
,
CB
11
,
CB
12
,
CB
13
,
CB
14
,
CB
15
,
CB
16
,
CB
17
,
CB
18
,
CB
19
,
CB
20
,
CB
21
,
CB
22
,
CB
23
,
CB
24
,
CB
25
,
CB
26
,
CB
27
,
CB
28
,
CB
29
,
CB
30
,
CB
31
,
CB
32
,
CB
33
,
CB
34
,
CB
35
,
CB
36
,
CB
37
,
CB
38
,
CB
39
,
CB
40
,
CB
41
,
CB
42
,
CB
43
,
CB
44
,
CB
45
,
CB
46
,
CB
47
,
CB
48
,
CB
49
,
CB
50
}
Parameters: {KA,KB}

Defining the model input for NonlinearModelFit command
model[KA_,KB_][i_,t_]:=Through[pfunall[KA,KB][t],List][[i]]/;And@@NumericQ/@{KA,KB,i,t};
In[]:=
datanoisy=Flatten[Join[{Table[{50,datareal1[[i]][[1]],datareal1[[i]][[2]]},{i,41,Length[datareal1]}],Table[{100,datareal2[[i]][[1]],datareal2[[i]][[2]]},{i,41,Length[datareal2]}]}],1];
In[]:=
Length@datanoisy
Out[]=
242
In[]:=
fit=NonlinearModelFit[datanoisy,{model[KA,KB][i,t],1.0>KA>0.8&&1.4>KB>1.2},{KA,KB},{i,t},Method->"InteriorPoint"](*callingtheNonlinearModelFitcommandandgatheringrelevantstatisticaldata*)
Out[]=
FittedModel
model[0.98,1.38][i,t]

In[]:=
nlm=Normal@fit
Out[]=
model[0.979998,1.38][i,t]
In[]:=
Total[fit["FitResiduals"]^2]
Out[]=
7.77464×
-7
10
In[]:=
p3=Plot[{model[0.979998,1.38][50,t],model[0.979998,1.38][100,t]},{t,0,1600},PlotStyle->Gray];
In[]:=
Show[p1,p3]
Out[]=
CA
50
CB
50

CITE THIS NOTEBOOK

Recovering equilibrium constants from simulated noisy chromatography data​
by Housam Binous​
Wolfram Community, STAFF PICKS, August 30, 2026
https://community.wolfram.com/groups/-/m/t/3789514