This post contains an implementation of the IWA's Activated Sludge Model No. 1 (ASM1) in octave.
Models are a widely-used tool in the design and operation of wastewater treatment systems. For my job, I use a commercial simulator, but I thought it would be a useful exercise to write an implementation of a simpler model (i.e. ASM1) on my own—to better understand how these models work "under-the-hood". Most wastewater models have a similar structure, built around the "Gujer matrix", which relates the rates of various processes to the changes they cause on different components.
Because the model was already defined and expressed in matrix form, I decided to code a script for it in octave. Octave is an open-source computer language that is highly similar to Matlab; vectors and matrices are its default data type, so it is well-suited for this application.
My script is included below. It is written for batch kinetics, so the mass balance only includes transformations of the initial components, not any flow in or out of the reactor. The dissolved oxygen concentration is fixed (as in a respirometer or a system with DO control) at 2.0 mg/L. The value could be changed to zero to simulate denitrification.
Following the code, I've included an example of its use with input values and a few output graphs.
Reminder: This blog is solely my own and may not reflect the views of my employer or association(s). If you use this code, you do so at your own risk.
The Code
% ASM1 (IWA's "Activated Sludge Model No. 1") has a stoichiometric matrix, S (8x13),
% a vector of process rates, Rho (8x1), and a vector of component conversion
% rates, r (1x13).
% S will be fixed in each run of the model while Rho and r depend on the current
% concentrations and will thus change at each time-step (dt).
%
%%% Calculation Procedure
% 1. Read in the initial concentrations, Ci (1x13 vector)
% 2. Set the values of the various coefficients (i_, f_, Y_, k_, etc.)
% 3. Initialize the stoichiometric matrix, S
% 4. For an appropriate time-step, dt, do the following calculations
% * Use Ccurr to determine Rho
% * r = Rho' x S; (verify with rules of matrix multiplication???)
% * Ccurr = Ccurr + r x dt;
% * For Oxygen, Ccurr(8), it may be better to modify to have continuous input and record the uptake rate
% * At a suitable interval (Interval), record Ccurr in Ct (N x 13)
% 5. When the simulation is finished, Ct can be used for output. A plot should be produced.
%
%% NOTE: This script is only set up for batch kinetics in a single reactor.
%
% Table numbers refer to IWA Scientific and Technical Report No. 9.
%
%% Read-in initial concentrations
Ci = zeros(1,13);
fprintf("Enter intial concentrations of each model component, in mg/L.\n");
Ci(1,1) = input("Soluble inert organic matter, Si: \n");
Ci(1,2) = input("Readily biodegradable substrate, Ss: \n");
Ci(1,3) = input("Particulate inert organic matter, Xi: \n");
Ci(1,4) = input("Slowly biodegradable substrate, Xs: \n");
Ci(1,5) = input("Active heterotrophic biomass, Xb,h: \n");
Ci(1,6) = input("Active autotrophic biomass, Xb,a: \n");
Ci(1,7) = input("Particulate decay products, Xp: \n");
Ci(1,8) = 2.0; % Assume a fixed O2 concentration and calculate the OUR.
Ci(1,9) = input("Nitrate- and nitrite- nitrogen, Sno: \n");
Ci(1,10) = input("Ammonia-nitrogen, Snh: \n");
Ci(1,11) = input("Soluble biodegradable organic nitrogen, Snd: \n");
Ci(1,12) = input("Particulate biodegradable organic nitrogen, Xnd: \n");
Ci(1,13) = input("Alkalinity - molar units, Salk: \n");
Ccurr = Ci;
%% Set the values of the coefficients and constants
% Default values from Table 5, p. 24
YA = 0.24; % g cell COD formed/ g N oxidized
YH = 0.67; % g cell COD formed/ g COD oxidized
fp = 0.08; % dimensionless, fraction of biomass going to inert products
iXB = 0.086; % g N in biomass / g COD in biomass
iXE = 0.06; % g N in endogenous mass/ g COD in endogenous mass
MuH = 6.0; % per day; heterotroph growth
Ks = 20.0; % g COD/m^3; COD half-saturation coefficient
KOH = 0.20; % g O2/m^3; oxygen switch
KNO = 0.50; % g NO3-N/m^3; nitrate switch
bH = 0.62; % per day; heterotroph decay
bA = 0.10; % per day; autotroph decay
Etag = 0.8; % anoxic growth discount
Etah = 0.4; % anoxic hydrolysis discount
kh = 3.0; % g slowly biodegradable COD/g cell COD/day; hydrolysis rate
KX = 0.03; % g slowly biodegradable COD/g cell COD; hydrolysis half-saturation coeff.
MuA = 0.80; % per day; autotroph growth
KNH = 1.0; % g NH3-N/m^3; ammonia half-saturation coefficient
KOA = 0.4; % g O2/m^3; oxygen switch
ka = 0.08; % m^3/g COD/d; ammonification rate
%% Initialize the stoichiometric matrix
S = zeros(8,13);
%process 1: aerobic growth of heterotrophs
S(1,2) = -1/YH;
S(1,5) = 1;
S(1,8) = -(1-YH)/YH;
S(1,10) = -iXB;
S(1,13) = -iXB/14;
%process 2: anoxic growth of heterotrophs
S(2,2) = -1/YH;
S(2,5) = 1;
S(2,9) = -(1-YH)/(2.86*YH);
S(2,10) = -iXB;
S(2,13) = ((1-YH)/(14*2.86*YH)) - iXB/14;
%process 3: aerobic growth of autotrophs
S(3,6) = 1;
S(3,8) = -(4.57-YA)/YA;
S(3,9) = 1/YA;
S(3,10) = -iXB - (1/YA);
S(3,13) = -(iXB/14) - (1/(7*YA));
%process 4: decay of heterotrophs
S(4,4) = 1-fp;
S(4,5) = -1;
S(4,7) = fp;
S(4,12) = iXB - fp*iXE;
%process 5: decay of autotrophs
S(5,4) = 1-fp;
S(5,6) = -1;
S(5,7) = fp;
S(5,12) = iXB - fp*iXE;
%process 6: ammonification
S(6,10) = 1;
S(6,11) = -1;
S(6,13) = 1/14;
%process 7: hydrolysis of entrapped organics
S(7,2) = 1;
S(7,4) = -1;
%process 8: hydrolysis of entrapped nitrogen
S(8,11) = 1;
S(8,12) = -1;
%% Simulate with suitable time-step;
dt = 1/1440; %in days
Interval = 1/24; %in days
maxtime = input("How long should the simulation run (in days)? \n");
Record = 1; %index of Ct
Rho = zeros(8,1); %Create Rho and r as empty vectors
r = zeros(1,13);
for time = 0:dt:maxtime
%Record Ct
if mod(time,Interval)==0
Ct(Record,:) = Ccurr; %Save the concentrations in a matrix for output
Ct(Record,8) = r(1,8); %Replace the O2 concentration (constant) with the OUR in mg/L*d.
Timet(Record,1) = time; %Save the times in a vector for output
Record = Record + 1; %increment the index
end
%Calculate Rho
Rho(1,1) = MuH*(Ccurr(1,2)/(Ks+Ccurr(1,2)))*(Ccurr(1,8)/(KOH+Ccurr(1,8)))*Ccurr(1,5);
Rho(2,1) = MuH*(Ccurr(1,2)/(Ks+Ccurr(1,2)))*(KOH/(KOH+Ccurr(1,8)))*(Ccurr(1,9)/(KNO+Ccurr(1,9)))*Etag*Ccurr(1,5);
Rho(3,1) = MuA*(Ccurr(1,10)/(KNH+Ccurr(1,10)))*(Ccurr(1,8)/(KOH+Ccurr(1,8)))*Ccurr(1,6);
Rho(4,1) = bH*Ccurr(1,5);
Rho(5,1) = bA*Ccurr(1,6);
Rho(6,1) = ka*Ccurr(1,11)*Ccurr(1,5);
Rho(7,1) = kh*((Ccurr(1,4)/Ccurr(1,5))/(KX+(Ccurr(1,4)/Ccurr(1,5))))*((Ccurr(1,8)/(KOH+Ccurr(1,8)))+Etah*(KOH/(KOH+Ccurr(1,8)))*(Ccurr(1,9)/(KNO+Ccurr(1,9))))*Ccurr(1,5);
Rho(8,1) = Rho(7,1)*Ccurr(1,12)/Ccurr(1,4);
%Calculate r and Ccurr
r = Rho' * S;
Ccurr = Ccurr + r*dt;
Ccurr(1,8) = Ci(1,8); % reset O2 concentration
end
%% Plot the results
tCOD = Ct(:,1) + Ct(:,2) + Ct(:,3) + Ct(:,4);
MLVSS = (Ct(:,3) + Ct(:,4) + Ct(:,5) + Ct(:,6) + Ct(:,7))/1.42;
OUR = -1*Ct(:,8);
plot(Timet,tCOD,'r-');
title('Total COD');
figure;
plot(Timet,MLVSS,'g-');
title('Biomass');
figure;
plot(Timet,OUR,'b-');
title('Oxygen Uptake Rate');
%Type 'figure(1)', 'figure(2)', 'figure(3)' to cycle back through figures
%With a figure selected, save it using the command 'print('filename.png','-dpng')'
Input and Output
Here are some graphs produced from the script above, together with the input values I used:
(Units are days and mg/L)
On a different topic...
In the simplified mark-up language that ghost (i.e. this blogging platform) uses, blocks of code should be preceded by 4 spaces. These can be added efficiently in a bash terminal using sed: cat ASM1_batch2.m | sed 's/^/ /'