Auto-correlation robustness of factorial designs and GAMS in studying the effects of process variables in a dual-objective adsorption system

The performance of factorial designs is still limited due to some uncertainties that usually intensify process complexities, hence, the need for inter-platform auto-correlation analyses. In this study, the auto-correlation capabilities of factorial designs and General Algebraic Modeling System (GAMS) on the effects of some pertinent operating variables in wastewater treatment were compared. Individual and combined models were implemented in GAMS and solved with the trio of BARON, CPLEX and IPOPT solvers. It is revealed that adsorbent dosage had the highest effect on the process. It contributed the most effect toward obtaining the minimum silica and TDS contents of 13 mg/L and 814 mg/L, and 13.6 mg/L and 815 mg/L from factorial design and GAMS platforms, respectively. This indicates a concurrence between the results from the two platforms with percentage errors of 4.4% and 0.2% for silica and TDS, respectively. The effects of the mixing speed and contact time are negligible.


Introduction
In process systems engineering, uncertainty, an inexorable concomitant of unknown events and imprecise information exists in many processes, and this usually influences the decision-making process and system sustainability, as well as intensifies the process complexities. Examples of the uncertainties that exist in every process designs and operations are: platform uncertainty, which usually reflects the identification of the best fitting design for a particular process, and the principles upon which the different platforms utilized for process designs are built; and process uncertainty, reflecting the uncertainty faced in determining the requirements that the control variables must fulfill in arriving at the best process besides their interactive effects on the objective functions. Once the main sources of a process uncertainty having significant effects are identified, then only those variables will be used in the final process design. In the design of experiments intended for bigger and complex systems, the larger number of involved variables and experimental runs makes certain analyses impractical or prohibitive. For any process to be sustainable, the development of sound approaches that characterize the uncertainty within process variables as well as provide decision-makers with a comprehensive analysis of the effects of variable variations on the process responses is crucial. To tackle such uncertainty issues, several mathematical modelling ( One of the broadly applied statistical methods in design of experiments (DoE) is the factorial design which is usually employed in studying the interactive effects of two or more variables/factors on an objective function or response variable. The full-factorial (FF) design is an authoritative procedure of studying the effects of numerous control variables with multiple levels on a response/objective function. For instance, if the effect of a control variable on the response depends on the level of another control variable(s), then there is an interaction between them. Thus, through factorial designs, we can study not only the core effects but also the interactive effects between variables. Furthermore, the technique involves the statistical experimental design which is observed as an efficient way of obtaining the exhaustive amount of information with the least sum of experiments probable (Zhang et al. 2009).
The objective/response of the FF design is anticipated to be seamlessly modeled by a function (linear) of the control variables, and then, a first-order model will be the approximation function. Moreover, since there are typically just two levels for each variable, the usual assumption is to have an approximately linear objective function over the selected levels of variable range. This especially holds for the study of a system/process just starting with a necessity to conduct variable screening experimentation. The utility of factorial designs is imperative since they are extensively employed in research works to form the basis for the decision frameworks or designs of highly substantial practical value. The regression model for the prediction of an objective function (Y) for such designs is given in Eq. 1 (Adeniyi et al. 2019a;Montgomery 2004): The regression model(s) generated for the objective function(s) besides ANOVA tests are employed in the decisive evaluation of the statistically significant variables or interaction (Adeniyi et al. 2019b). However, the performance of factorial among other methods is still limited due to some uncertainties that usually intensify process complexities of which conflicting regions of variable interest are a basic example (Amosa and Majozi 2016).
General Algebraic Modeling System (GAMS) is an optimization and modeling framework whose concepts and designs are drawn from both interactive database theory and mathematical programming. It consists of a language compiler and integrated highly efficient solvers for formulating, solving, and/or analyzing an optimization problem. GAMS could well be adapted for complex, large-scale modeling applications, and can aid in building sustainably huge models that may rapidly be well adapted to newer situations. The platform is designed for modeling linear, nonlinear and mixed integer optimization problems, and more interestingly, complex and large types of optimization problems requiring many revisions to establish an accurate model can be handled (GAMS 2015). It has been successfully applied in a variety of applications such as optimizing wastewater treatment networks (Galan and Grossmann 1998), membrane networks synthesis and optimization (Abass and Majozi 2016;Buabeng-Baidoo and Majozi 2015;Khor et al. 2011;Mafukidze and Majozi 2016) and water network synthesis as excellently reviewed by Khor et al. (2014). Besides, GAMS has been successfully employed for multi-objective models of the stacks of thermo-acoustic engines/refrigerators (Tartibu et al. 2014(Tartibu et al. , 2015. Recently, the removal of silica, TDS and other pollutants from industrial wastewater has become a thought-provoking task for the environmental process engineers Ighalo and Adeniyi 2020b). High silica and TDS contents of wastewater will always be a concern, especially when the wastewater is to be reclaimed for specific reuse purposes. Among other possible damages to reactors, these two constituents can result in corrosion, scale deposition, fouling in membranes and boiler water solutes in steam if they are not kept within the necessary threshold range Jami et al. 2013). Several processes have been adopted for removal of silica and TDS (Basha et al. 2008;Broséus et al. 2009;Latour et al. 2014Latour et al. , 2015Lee et al. 2008;Sik Ali et al. 2004;Weeks et al. 1992). Overall, some of these methods are efficient, but each one of them has their limitations ranging from process complexities to cost implications. This creates the room for a need (1) Y= 0 + 1 A + 2 B + 3 C + 12 AB + … + to develop innovative and cheaper technologies that can remove TDS at a larger scale with a focus on minimizing energy and material costs, maximizing water recovery, without any compromise on the environment (Ighalo and Adeniyi 2020a). Having realized that there is a dearth of documentation about how efficient adsorbents could be in removing silica and TDS, this study has focused on adsorption process using a locally produced adsorbent derived from a waste called oil palm empty fruit brunch. Specifically, there is a need to investigate the interactive trends between the control variables and the pollutants compared to other contaminants.
In our previous investigation, we established how regression models developed from response surface methodology (RSM) could be well supported and facilitated by GAMS, especially when the regions of variable interests are conflicting (Amosa and Majozi 2016). Our main aim and contribution in the present study were to evaluate two main strategies, i.e., the comprehensive study of the variable interactions and/or the effects of process variables on the uptake of silica and TDS from high strength wastewater using FF design. The second evaluated strategy, which to the best of our knowledge has not been previously reported, is to autocorrelate the variable effects of the two platforms (factorial design and GAMS) toward making a steadier trade-off decision on the appropriate choice of the variable(s) having the most significant effect on the objective function(s). This will add to the knowledge necessary for choosing the best adsorption operating conditions towards achieving local optimum, and economical utility of adsorbents and energy consumption concerning the mixing requirements for silica and TDS removal from wastewaters of similar chemistry. It was also expected that this work would give a justifiable basis for the accurate evaluation of interactive effects between the specific variables and objective functions in similar adsorption processes. Also, it will form a basis on which large adsorption systems can be designed and controlled using minimal costs, time and space. Discussion on adsorbent characterization, and kinetics and equilibrium studies is, however, outside the scope of the current study.

Adsorbent
A locally produced low-cost adsorbent (powdered activated carbon-PAC) with steam pyrolysis activation route was utilized for the uptake of silica and TDS. The comprehensive features of the adsorbent are hereby presented in Table 1 and have been exclusively reported elsewhere (Amosa 2015;Amosa et al. 2016bAmosa et al. , 2014Amosa et al. , 2015Amosa et al. , 2016c.

Wastewater sampling and preservation
The wastewater (biotreated POME) employed in this study was sampled from the Sime Darby Mill located in Carey Island, Malaysia. The physicochemical characterization of the sample was carried out using the Standard Methods for the Examination of Water and Wastewater approved by the American Public Health Association (APHA) and HACH Company (APHA 2005; HACH 2012). Besides, the sample was appropriately preserved strictly implementing the US EPA preservation standards (Alam et al. 2016;US EPA 1982). The profile of the measured contaminants in the wastewater is presented in Table 2.

Water quality characterization
Silica content was determined with the aid of a HACH spectrophotometer model DR 5000 (HACH, Loveland, USA), while TDS content was determined with HACH sensION™ 7 model (HACH, Loveland, USA). Extensive procedural methods for the measurement of all contaminants have been reported elsewhere (Amosa 2015;Amosa et al. 2014Amosa et al. , 2015Amosa et al. , 2016c.

Software packages
The design of experiment (DoE) for the process execution and statistical analyses was implemented in Design Expert® software version 10.0.4.0 (Stat-Ease Inc., MN, USA). The first-order polynomial models that were generated were directly implemented in GAMS without any adjustment to permit easy comparison of the results. Meanwhile, the models are nonlinear; hence, they are handled as NLP models implemented in GAMS version 24.4.3 (GAMS 2015), and solved with BARON (Tawarmalani and Sahinidis 2005), together with CPLEX and IPOPT solvers, for comparison. The three different solvers were applied to see if there is any discrepancy in the results obtained. A personal computer characterized by a Core i5-3210 M CPU processor at 2.50 GHz and 12.0 GB RAM was utilized to carry out all computations. The three control variables were bounded as presented in Eq. 2, representing the necessary constraints for the solver in searching for the best solutions in the regions set for the constraints. This set of constraints is utilized for all models in this study.
Since real wastewater is employed in this study, and the main objective is to monitor the removal of the contaminants contained therein, the final concentrations of silica and TDS are set as the objective functions. The variables A (adsorbent dosage), B (mixing speed) and C (contact time) in the constraints are employed to identify the vital variables that affect the adsorption process. This enhances the discovery of the ideal combination of categorical variables. All the results were compared using the percentage error approach. Percent error, which is the absolute value of the difference of two values divided by the theoretical value, was utilized here because the experimental result was compared with a theoretical (model) value.

Batch adsorption tests using 2 3 full-factorial design (first-order model)
To observe the manners in which the control variables and their interactions are influential to the objective functions, a full-factorial design was performed using the Design Expert® software. Batch adsorption experimentation was conducted to evaluate the adsorption efficiency and the effects of operating variables using a 2-level FF design. To be certain of the homogeneity of the representative sample throughout the experimentation, the sample was homogenized in large aliquot to abate random sampling errors due to segregation effects. The agitated system is used to maintain balanced quantities of the contaminants in different phases based on concentration levels. And since soluble solids are involved, the agitation is used to increase interaction between the particles and avoid uneven accumulation at one point. Though turbulent regions of the homogenized solution may be chaotic and difficult to determine, the agitation process enhanced the wastewater samples to be taken from well-mixed regions.
The experimentations were carried out by mixing varying doses of the adsorbent with 100 mL of the wastewater sample, followed by agitation in 250-mL-sized Erlenmeyer flasks at various mixing speeds and contact time, depending on the set of variable levels in the factorial design. The temperature and pH values of the sample were left unchanged to specifically simulate the nature of the ready-to-discharge wastewater. The wastewater profile is characterized using the standard methods (APHA 2005; HACH 2012) before and after the adsorption processes. All the experiments are triplicated, and the average values are logged. Throughout the adsorption process in this investigation, the independent/control variables: adsorbent dosage, mixing speed and contact time are, respectively, represented as A, B and C. The levels and range of the control variables involved in the experimental design are shown in Table 3. The values for the variable range were adjusted based on the results of preliminary study and previous findings (AbdulHalim et al. 2011;Alam et al. 2016Alam et al. , 2009Alkhatib et al. 2011;Amosa et al. 2014;Emad 2010;Neşe and Ennil 2008;Razali et al. 2010;Tumin et al. 2008).

Analyses of the effects of control variables on responses
Based on the 2-level FF design, monitoring of the adsorption trends of each constituent was performed to obtain the variable relational effects on the adsorption process through a combination of 11 experimental runs as offered by the design software. The experimental design and results for the removal of silica and TDS are presented in Table 4. The results are appropriately analyzed through the ANOVA. ANOVA implies a statistical algorithm that splits the overall   Turan and Ozgonenel 2013). In this study, the statistical significance of the variation in mean square ratio due to the residual errors and regression of the mean square are appraised with the application of ANOVA . It is observed that some of the control variables showed significant effects (i.e., p ˂ 0.05) while some of them are not significant (i.e., p ˃ 0.05)-see Tables 5 and 6.

Silica removal
With an initial concentration of silica (73 mg/L), the lowest residual concentration of 13 mg/L was obtained at the operating conditions of 5 g adsorbent dosage, 100 rpm mixing speed and 60 min of contact time. This indicates an 82.2% removal efficiency from the adsorption process. In Table 5, the ANOVA representing the model that adjusts for the significant curvature on silica removal is presented. The presented ANOVA is based on the default model used for diagnostic plots as well as the prediction and model plots. As observed in Table 5, it could be inferred from a model F-value of 453.38 that the model is significant, and there is a chance of 0.22% that "Model F-value" this big could occur because of noise. A p-value of 0.0022 for the overall model signifies the strong significance of the model, since the values of "Prob. > F" value of ˂ 0.050 indicate that the model terms are highly significant and values greater than 0.100 signify insignificant model terms. The "Curvature F-value" of 1230.01 suggests that in the design space, there is a significant curvature (measured by a difference between the center points average and factorial points average). Also, a chance of only 0.08% exists for a large "Curvature F-value" like this to occur as a result of noise. The signal-to-noise ratio is appropriately measured by adequate precision, and any ratio greater than 4 is deemed desirable. The signal-tonoise ratio as given by the model is 65.105 which implies an acceptable signal.
The relationship between the predicted and actual values is depicted in Fig. 1, and the linear trend proves that the two sets of values do relatively concur with each other. In this case, the significant model terms appear to be variables A, B, C, alongside the AC and ABC terms with R 2 and adjusted R 2 values of 0.9994 and 0.9972, respectively. Such magnitudes of R 2 values are clear indications that the model successfully correlates the objective functions to the variables studied. Furthermore, this is an indication that the model is fit to be utilized in the navigation of the design space. The model has an SD of 0.58, and the final regression model equation was generated by the Design Expert® software for the lowest concentration of silica residual in terms of the actual and coded variables as, respectively, given in Eqs. 3 and 4: Moreover, 3D variable interactive plots are also possible since the study conducted here involves at least two control variables. The 3D plots (Fig. 2) depict the interaction of the three variables that are involved in the silica sorption process. Figure 2a shows the relationship between the mixing speed and adsorbent dosage. An increase in the dosage of 2-5 g resulted in higher silica uptake from the effluent as evident from the low trend of silica residuals. Though not as highly significant as that of the dosage, it was also observed that as the mixing speed increases from 100 to 200 rpm, there is a slight reduction effect on the silica residuals. In the illustrated 3D diagram shown in Fig. 2b, the contact time and dosage had a significant effect on the process. The silica residual tends toward a rapid reduction as both variables increase, and the presence of curvature in the surface plot indicates a significant interaction between the variables. 3D plot in Fig. 2c shows that contact time and mixing speed interaction affected as the silica concentration continues to slightly decrease evident from the slight imbalances exhibited in the 3D plot axes.  It also helps to view the results in the cubical factor space. Figure 3 shows a cubical plot showing the combined effects of all the three control variables on the objective function. Since the variables of interest here are A, B and C, they were picked by default by the software program. The shown values are the predicted values, which permits plots to be made even with missing actual data. It is observed that silica is minimum at the A + , B − , C + settings (lower back right corner of the cube).

GAMS validation for minimum silica content possible
The model Eq. 3 was subjected to the constraints in Eq. 2 and solved with the BARON, CPLEX and IPOPT solvers in GAMS platform to find the local optimum for silica removal.
The solver utilizes the set upper and lower boundaries to determine the level and marginal values. The three solvers interestingly gave similar results of local optimum, levels and marginal values as shown in Table 6. Furthermore, the platform gave a model statistics of 11 non-zero elements, 6 single variables, a derivative pool of 20, 6 nonlinear N-Z, 29 constant pool, a code length of 52 with a generation time of 0.375. The level represents the proposed local optimum given as 13.6 mg/L but at the operating settings of 5 g of adsorbent dose, 100 rpm of mixing speed and 60 min of mixing time in which concurs with the results from the factorial design analysis.
Moreover, the interesting part of the analysis is given in the marginal column. This has a certain meaning that the set constraints are the binding constraints. It means that if the lower limits are changed, then the objective function value will change (Soroudi 2017). The marginal values show the sensitivity coefficients of objective functions to the constraints (bounded control variables or equations). The marginal value represents the amount and direction of change in the objective value given a unit increase in the active variable bound (GAMS 2018). This will enable the decision-maker to identify the binding constraint having the most influence on the objective function (Soroudi 2017).
The minimum value obtained in the GAMS platform represents a synergistic percent error of 4.4% when compared with the experiment-based result (13 mg/L). This solution can be physically interpreted thus: • Keeping the adsorbent dosage and the mixing speed at levels suggested by both the experimental and theoretical solutions ensures that the objective is adequately minimized. • Still, there exists a synergistic effect of the control variables on the response since the minimum values obtained from each platform are very close with a percent error of just 4.4%.

Total dissolved solids (TDS)
The minimum TDS residual of 814 mg/L (from an initial concentration of 1207 mg/L) was observed at the operating conditions of 5 g adsorbent dosage, mixing speed of 100 rpm and 60 min of contact time. This leaves the highest removal efficiency at 32.6%. Besides, 829 mg/L and 833 mg/L were the closest TDS residuals observed but at different operating conditions. As presented in Table 7, the ANOVA gives a model F-value of 3168.64 for the adsorption process thereby signifying a relatively high significance of the model. Also, there exists just a chance of 0.03% that such a relatively large model F-value could occur from noise effect. The general model has a p-value of 0.0003, thereby indicating that the model is highly significant since "Prob. > F" values of < 0.05 are known to indicate significant model terms. The "Curvature F-value" of 55.68 infers that in the design space, there is certainly a significant curvature and the ANOVA suggests that for such a large curvature F-value to occur due to noise, there exists only a chance of 1.75%. An adequate precision ratio of 155.103 indicates an adequate signal since a ratio greater than 4 is deemed desirable. Figure 4 shows a plot depicting the correlation between the predicted (model) and actual (experimental) TDS residual concentration values. The linearity of the correlation plot proves that a relative agreement between both values exists. In this case, variables A, B, C and BC term represent the significant model terms with R 2 and adjusted R 2 values of 0.9999 and 0.9996, respectively. Such high R 2 values are indications that the model has successfully correlated the studied variables and the objective response. An SD of 1.15 is exhibited by the model which is an indication that in the design space, it can be used for navigation. Valid mathematical relationships of the variables and objective function are obtained for TDS residual concentration as a result of the ANOVA analyses of the experimental data selected via the FF design. The regression models generated by the Design Expert® software in the form of actual and coded variables are presented in Eqs. 5 and 6, respectively.
In this case as well, 3D plots (Fig. 5) depicting the interactive behaviors between the objective function and all control variables are generated by the software for analyses. The plots depict the effects of operating factors' interactions on the abatement of TDS. The interaction between the dosage and contact time is noticeable as the increase in contact time and dosage from 30 to 60 min and 2 to 5 g, respectively, favored TDS abatement. However, 3D plots in Figs. 5a and c suggest favorable TDS abatement trend when any of the factors interacted with decreasing mixing speed from 200 to 100 rpm as its increase did not favor TDS removal. This is not unprecedented as a similar observation on TDS sorption by activated carbon was also reported by Mortula and Shabani (2012). With the inspection of Fig. 5b, it is evident that any increase in the adsorbent dosage resulted in a lower concentration of TDS. This could be attributed to the fact that increased dosage provided extra surface area, available pores and functional groups ready for the TDS entrapment. And since the activated carbon is a porous material that acts as a sieve for the contaminants, it could be inferred that the PAC may have acted as a filter medium. It has been reported that adequate filter media such as porous materials including membranes, activated carbon, zeolitic materials or roughing filters are usually more efficient in TDS removal (AquaFit4Use 2010; Karakulski et al. 2006;Karimi et al. 2015;Nkwonta and Ochieng 2010;Pal et al. 2014;Song 2004;Wu et al. 2007). Besides, it must be noted that the variable effects on TDS also give the idea of the behaviors of salinity and electrical conductivity under similar operating conditions, since the three parameters are all closely related. An increase in any of the parameters leads to an increase in the other two parameters and vice versa.  Moreover, if the cube plot (Fig. 6) is also inspected, it is observable that the lowest TDS concentrations are achieved with A + , B − and C + which also corresponds to the same operating conditions attained in the case of silica concentration. Generally, it is obvious that the obtained model values in all cases were close to the experimental values and that is indicative of the fact that the developed models are efficacious in unfolding the correlation between the variables and the objective functions. Also, according to the ANOVA, the most significant of all the control variables is the dosage in both cases due to their very low p-values compared with the other two variables. This is indicative of the fact that a slight change in the dosage will significantly affect the two responses or objective functions. The model qualities in each case are also assessed based on the coefficient of determination (R 2 ) values with all R 2 values observed to be very close to unity. Besides, it is also detected that R 2 values are in reasonable agreements with their respective adjusted R 2 values for all both pollutants studied, and very low SD values are exhibited by the models. As R 2 values tend to unity with very low SD, the model becomes better indicating an established resilience of agreement between the predicted and actual values of the objective functions.

GAMS validation for minimum TDS content possible
Model Eq. 5, which is based on the actual variables and response values, was also subjected to the constraints in Eq. 2 and implemented in GAMS. The result proposed a local optimum setting of 5 g of adsorbent dose, 100 rpm of mixing speed and 60 min of mixing time as with a minimum TDS content of 815.5 mg/L. Though this value is a bit higher than the minimum value (814 mg/L) obtained from the experimental design, both values are still synergistic since the percentage error is only 0.18%. The result can thus be physically interpreted as follows: • Maintaining the levels of adsorbent dosage and mixing time as that obtained from both platforms is sufficient to obtain the minimum TDS content. • The control variables led to a synergistic effect on the objective as the residual concentration was 0.18% higher than the value provided by the model when GAMS was not implemented.

Combined solutions of the objective functions under factorial design and GAMS platforms
It was observed that the pollutants experienced their lowest residual concentrations (13 mg/L silica and 814 mg/L TDS) at the same non-conflicting operating conditions of 5 g of adsorbent dosage, 100 rpm of mixing speed and 60 min of contact time (Table 3). However, with more than two variables, the best or optimal conditions have to be obtained using numerical optimization methods as suggested by the software. This will allow the simultaneous consideration of all objective functions involved in the design and find a single operating setting and value of the objective function. The statistical analysis of the Design-Expert® software version 10.0.4.0 offered 28 numerical solutions for the best removal of silica and TDS with desirability values ranging from 0.747 to 1.000 for their possible minimum contents. The application of desirability function to systems in designs of experiments has been fully elucidated elsewhere (Amosa and Majozi 2016). Experimental validation of operating conditions having the highest desirability value of 1.000, and operating setting of 5 g of dosage, 100 rpm of mixing speed and 60 min of contact time was carried out, and it was observed that the validated results fully agreed with the numerically proposed results. Subsequently, to simultaneously find the minimum possible residual contents of the two pollutants, the process objective function models Y 1 and Y 2 , based on the actual values, were combined as a single dual-objective function (DOF) and represented as R for finding the best operating process conditions as follows: The DOF is subjected to the constraints in Eq. 2, and as GAMS approach was implemented in solving the DOF, a solution of 13.6 mg/L of silica and 815.5 mg/L of TDS was predicted with the set of control variables: 5 g of dosage, 100 rpm of mixing speed and 60 min of mixing/contact time.   Fig. 6 Cube plot of the interactions between the control variables and TDS uptake Moreover, the model statistics of the GAMS-based solution revealed that there are 3 single and blocks of equations, 6 single and blocks of variables of which 3 variables are projected, 11 nonzero elements, 6 nonlinear N-Z, a derivative pool of 20, constant pools of 29 and a code length of 52. More importantly and just as the ANOVA analyses suggested for both uptake cases, GAMS platform also suggests based on the marginal values that a slight change in the dosage will significantly affect the objective function by a marginal value of − 49.126 (see Table 6). Other two control variables (mixing speed and contact time) are not as significant as the dosage due to their very low and negligible marginal values. It appeared that both the GAMS-and factorial-based solutions agreed with the experimental values as the % errors were all less than 5%. Moreover, GAMS implementation in simultaneously solving the DOF problem appeared more economically viable since only a single set of control variables was suggested as the optimum. This was easier to validate experimentally even with triplicates, coupled with the fact that the optimal was somewhat close to that obtained from the factorial-based solution.

Conclusions
The statistical regression models for the two scenarios have been developed and their suitability confirmed with the ANOVA for prediction of the effects of variable interactions and interest regions for the process. Subsequently, the model is confirmed with three (3) solvers under the GAMS platform. Conclusions of the present investigation are abridged as follows: 1. The use of GAMS in combination with factorial design through ANOVA and marginal effects proved an effective and insightful way to characterize the relationships between control variables representing the constraints on the objective functions. 2. There is a close agreement between the results obtained from the two platforms regarding the variable significance, and both pointed to the fact that dosage is the most significant variable in the process. A slight change in the variable leads to a significant change in the objective functions. 3. This study further supports the use of the statistical design approach for the model development in terms of variable interactive effects that contributes to the suitability of the desired changes in the process being considered. 4. The studied interacting effects from the model are examined, and these expressed the close similitude with experimental interpretations. The results and regions of variable interests proposed by both the factorial and GAMS platforms concur with the experimentally derived results. 5. The model is confirmed under GAMS platform in terms of the best region of process variables necessary for the sustainability of the process. 6. Despite the multicomponent nature of the wastewater, the process efficiency was 82.2% for silica removal, while a lower efficiency was observed in the case of TDS removal. 7. Results of the variable interactive effects and probable levels form the basis for a proposed decision framework.