## 3 Results and Discussion  
### Simulation Setup  
The proposed mathematical model consists of a set of five ODEs and two algebraic equations (Equations (1)-(9)) thus represented by seven dependent state variables \(x=[X_{1}\ S_{1}\ A_{C1}\ Q_{H_{2}}\ X_{2}\ A_{C2}\ Q_{CH_{4}}]\) and nine unknown parameters divided into two groups: \(p_{1}=[\mu_{1 max}\ K_{S_{1}}\ Y_{1}\ Y_{2}\ Y_{H_{2}}]\) and \(p_{2}=[\mu_{2 max}\ K_{S_{2}}\ Y_{3}\ Y_{CH_{4}}]\). Each individual of the population (chromosome in GA, firefly in FA, nest in CS, and coyote in COA) represents a potential solution to the problem--the corresponding unknown model parameter vector \(p_{1}\) or \(p_{2}\). Based on [54,55] and the authors' expertise, each parameter of the model (1)-(9) has been coded in the specific range (lower bound (\(Lb\)) \(\leq\) parameter \(\leq\) upper bound (\(Ub\))) as follows:  
\[ \begin{array}{l}0.01\ \leq\ \mu_{1 max}\ \leq\ 0.8;\ 0.1\ \leq\ K_{S_{1}}\leq\ 50;\ 0.1\ \leq\ Y_{1}\ \leq\ 15;\\ 0.001\ \leq\ Y_{2}\ \leq\ 40;\ 0.01\ \leq\ Y_{H_{2}}\ \leq\ 15;\\ 0.001\ \leq\ \mu_{2 max}\ \leq\ 0.8;\ 0.0001\ \leq\ K_{S_{2}}\leq\ 5;\ 0.1\ \leq\ Y_{3}\leq\ 50;\ 0.1\ \leq\ Y_{CH_{4}}\leq\ 10. \end{array}  \tag{10}\]  
Each initial solution for all four metaheuristic algorithms has been generated based on  
\[x_{j}=rand*(Ub-Lb)+Lb. \tag{11}\]  
Further improvement of the solution has been made based on the specific algorithm inspired by the considered behaviour in nature in the case of GA, FA, CS, and COA. New solutions are generated taking into account the bounds (constraints) of the model parameters. In addition, each newly improved solution is tested to see if it is within the defined bounds. The metaheuristic algorithms have been executed on Intel(r) Core(tm)7-8700 CPU @ 3.20 GHz, 3192 MHz (ASBIS Bulgaria Ltd., Sofia, Bulgaria), 32 GB Memory (RAM), Windows 10 pro (64 bit) operating system. Matlab and Simulink R2019a environment are used. Based on the presented pseudocodes (Algorithms 1-4), the algorithms are implemented in Matlab code. Mathematical model of corn steep liquor TSAD process (Equations (1)-(9)) is presented as a Simulink file (see Figure 1). Solver options are fixed-step size of 0.01 and ode4 (Kunge-Kutta) with TIMESPAN = [0 12]. To achieve the best metaheuristic algorithm performance, necessary adjustments of the GA, FA, CS and COA parameters depending on the problem domain have been provided [23]. Starting from an already explored range of parameters variation in the case of GA [29; 56; 57], FA [19; 34; 35; 36; 38], CS [21; 42; 43; 50], and COA [44; 45; 46] and taking into account the characteristics of the particular identification problem here, several pre-tests (based on trial-and-error method) have been performed to tune the parameters of GA, FA, CA, and COA. The four algorithms have been run for 100 iterations, which preliminary tests showed was sufficient for the algorithms to converge. The main GA functions are as follows: fitness function--linear ranking, selection function--roulette wheel selection, crossover function--double point crossover, and mutation function--bit inversion. A binary chromosome representation with a precision of binary representation 20 is used. The set of optimal values used for the main parameters of the algorithms after a series of adjustments is listed in Table 1.  
To ensure a fair comparison of the considered stochastic metaheuristics, each algorithm is run 30 times for 100 iterations; the initial solution is generated based on Equation (11) and the model parameter constraints (_Lb_ and _Ub_) are as presented in (10). A two-stage identification scheme is used.