Objective functions to maximise the recovered biogas, to maximise the
recovered N and P nutrients within digestate, and to minimise the aeration
requirement as a proxy to total operation cost were optimised individually to
observe the trade-offs. These objective functions were of the form shown in
Equation 34 for $k=6,7,8$. The objective function to maximise nutrient
recovery modelled the quality of the recovered digestate in terms of the total
nutrient (N and P) content in kg t-1. The biogas recovery objective function
modelled the total methane recovered, in m3 d-1, from both the UASB and the
sludge digester. Finally, the aeration requirement was modelled as the total
oxygen flow, in m3 d-1, required to reach the desired dissolved oxygen set
point in the aerobic zones.

$\min\sum_{k}\gamma_{k}\left(\hat{\bar{y}}_{k,v}\sigma_{v}+\mu_{v}\right)$
(34)

The $\epsilon$-constraint method was also adopted to obtain additional non-
dominated Pareto optimal solutions (Equation 35) where $\epsilon_{v}$
constraints output variable $v$ to be at least as good as a given upper bound.

$\sum_{k}\gamma_{k}\left(\hat{\bar{y}}_{k,v}\sigma_{v}+\mu_{v}\right)\leq\epsilon_{v}$
(35)

Solutions to the MINLP problem thereby highlight optimal process
configurations (optimised binary variables) as well as optimal process designs
(optimised continuous variables). Initially, the influent COD was treated as
an optimisation decision variable so as to demonstrate the multiple decision
criteria trade-offs. Subsequently, the influent COD was treated as an
uncertain parameter within a stochastic programming framework. All of the
MINLP problems were solved using BARON version 22.9.30 [111] with a maximum
time limit of 500 s using a CPU with 2.8 GHz Quad-Core Intel Core i7 processor
and 16 GB of RAM.

####  3.2.6 Stochastic programming formulation

To optimise the superstructure configuration under the uncertainty in the
influent COD concentration, the constraints and objective functions were
reformulated to express the expected values of the relevant surrogate model
outputs. For example the effluent quality constraints were reformulated to
constrain the expected COD, TN, and TP concentrations below the specified
quality limits whilst the expected biogas recovery, expected nutrient quality
of the digestate, and expected aeration requirement were specified as the
stochastic programming objective functions.

To this end, the influent COD was first redefined as an indexed parameter as
opposed to an optimisation decision variable where the probability of each
influent COD composition, $R_{s}$, being realised, $p_{s}$, was evaluated.
Practically, the realised influent COD values were taken as a subset from the
standardised sampled data $\bar{X}_{2}$. This simultaneously enables good
sample coverage of the probability distribution and realised parameter values.
The relevant constraints and objective functions were then expanded as a
weighted sum-product between the probability of each realisation, $p_{s}$, and
the realised surrogate model prediction of dependent variables,
$\hat{\bar{y}}_{k,v,s}$. Equations 36, 37, and 37 show the stochastic
programming reformulations of the effluent quality constraints for the
expected COD, TN, and TP concentrations, respectively.

$\frac{1}{\sum
p_{s}}\sum_{k}\gamma_{k}\left(\sum_{s}p_{s}\left(\hat{\bar{y}}_{k,3,s}\sigma_{3}+\mu_{3}\right)\right)\leq
50$ (36) $\frac{1}{\sum
p_{s}}\sum_{k}\gamma_{k}\left(\sum_{s}p_{s}\left(\hat{\bar{y}}_{k,4,s}\sigma_{4}+\mu_{4}\right)\right)\leq
10$ (37) $\frac{1}{\sum
p_{s}}\sum_{k}\gamma_{k}\left(\sum_{s}p_{s}\left(\hat{\bar{y}}_{k,5,s}\sigma_{5}+\mu_{5}\right)\right)\leq
5$ (38)

The surrogate models are thereby evaluated at $x_{2}=R_{s}$ where $x_{2}$ is
the input variable representing the influent COD (Table 11) and $R_{s}$ is the
influent COD value realised for scenario $s$ equivalent to the standardised
sampled values $\bar{X}_{s,2}$. The GPR constraint reformations are hereby
shown by Equation 39.

$\hat{\bar{y}}_{k,v,s}={}\left.\sum\limits_{i\in
n}\alpha_{i}\sigma_{f}^{2}\exp\left(-\frac{1}{2l^{2}}\left(\left(\bar{x}_{1}-\bar{X}_{i,1}\right)^{2}+\left(\bar{X}_{s,2}-\bar{X}_{i,2}\right)^{2}\right)\right)\right|_{k,v}$
(39)

Similarly, the classification NNs were reformulated such that the second input
dimension was parameterised using the sampled data as discrete scenario
realisations to yield predictions denoted as $\hat{y}_{k,9,s}$. Whilst, the
“choose exactly one” constraint for the discrete BNR pathways remained
unchanged (Equation 24), the classification feasibility constraint was updated
to ensure that solutions were robust to the uncertainties by ensuring the
process was feasible for all influent COD scenarios (Equation 40).

$\sum_{k}\gamma_{k}\hat{y}_{k,9,s}\geq 0$ (40)

Finally, the objective functions were updated to minimise the expected values
over the possible uncertainties as shown by Equation 41 for $v=6,7,8$.

$\min\frac{1}{\sum
p_{s}}\sum_{k}\gamma_{k}\left(\sum_{s}p_{s}\left(\hat{\bar{y}}_{k,v,s}\sigma_{v}+\mu_{v}\right)\right)$
(41)

##  4 Results and discussion

###  4.1 Computer experiment results

Figure 11 shows the samples, which were generated from the Python-GPS-X
interface in 374 s, distributed over the input space and the convergence
results. It can be seen that 7 out of the total 128 samples failed to converge
(about 5 %), primarily located at low UASB volumes at high COD loading for the
A2O, Bardenpho, and Johannesburg BNR pathways. The UCT process simulations
converged for all 32 input samples. Overall, these results demonstrate the
good sample space coverage exhibited by the static Sobol sampling over the 2
input dimensions. The total sampling time of 374 s equates to just less than 3
s per sample. Samples that failed to converge took roughly 3 times as long at
just under 9 s per sample which corresponds to the 3 allowed convergence
attempts per sample. Although the total sampling time here is acceptable, for
higher input dimensions necessitating greater numbers of input samples, and
for more complex simulations with greater convergence failure rates, the total
sampling time can contribute significant time costs to the surrogate modelling
framework. The efficient sample generation exhibited here coupled with the
high convergence rate indicates a well posed design space imposed by the input
variable selection and their corresponding bounds.

Figure 11: Computer experiments results showing well-spaced Sobol samples
within the input space. Converged samples are shown as black whilst non-
converged samples are shown as red. Additionally, training samples are shown
as filled whilst testing samples are shown as unfilled.

Figure 11 also shows the random split between training and testing data where
25% of samples were reserved for testing, resulting in 24 training samples and
8 testing samples. The same train-test split was applied to each discrete BNR
pathway so that the uncertainties introduced by this split were equivalent for
each discrete scenario. Specifically, Figure 11 shows how the random train-
test split interferes with the spacing of the 32 samples resulting in regions
of sparsely sampled training data within the 2-dimensional continuous search
spaces. Additionally, the resulting testing set is not well spaced which can
result in unrepresentative validation metrics for surrogate models trained to
represent the entire search space. The implications of the random train-test
split are highlighted in the subsequent sections.

###  4.2 Surrogate modelling results

####  4.2.1 Validation results

24 GPR models, each trained on 24 training samples in under 0.03 s, were
validated on the reserved testing data. The average MAEs across the 4 BNR
pathways as well as the mean output values are shown in Table 12. Whilst the
MAEs require knowledge of the magnitude of the underlying data, the MAPEs are
normalised such that a value of 1 represents an error of the order of
magnitude equivalent to the underlying data. In this way MAPEs can be directly
compared across variables of different magnitude as shown in Figure 12. It can
be interpreted that the GPR errors for effluent COD and TP concentrations are
high and that all other GPR errors are acceptable. It should be noted that
high errors are expected due to the random train-test split negatively
impacting the well-spaced properties of the training samples, and the
resulting placement of the testing samples in the most sparsely populated
regions.

Variable | Unit | MAE | $\mu$
---|---|---|---
effluent COD | g m-3 | 29.8 | 72.4
effluent TN | g m-3 | 2.7 | 9.6
effluent TP | g m-3 | 2.4 | 6.7
nutrient quality | kg t-1 | 2.5 | 36.3
biogas recovery | m3 d-1 | 67.8 | 707
aeration requirement |  $,000$ m3 d-1 | 7.5 | 38.5
Table 12: Average mean absolute errors and mean values for output variable
predictions across the 4 different BNR pathways.

The calculation of the MAPE also results in large values when the underlying
data values are small. For example, effluent concentrations of COD, TN, and TP
in the treated wastewater are very small for large part of the search space
(depicted later in this chapter in Figure 13). Additionally, the mean values
are often skewed by high effluent concentrations at poor performing process
design samples. The GPR validation results should therefore be used primarily
for greater interpretation of the optimisation results in conjunction with the
GPR predictive uncertainties.

Figure 12: GPR mean absolute percentage errors.

####  4.2.2 Response surface results

The response surfaces generated by the GPR models predicting effluent
concentrations of COD, TN, and TP, in addition to the feasible region
predicted by the classification NN trained on convergence labels, are shown in
Figure 13. The effluent quality response surfaces also depict the feasible
region separated by the contour placed at the maximum effluent concentration
limits. Together, the agglomeration of these 4 distinct feasible regions
defined the feasible search space for the subsequent superstructure
optimisation.

Figure 13 shows that the effluent COD quality limit results in infeasible
process designs at low UASB volumes for high influent COD concentrations.
Additionally, it can be observed that about 75% of the search space predicts
effluent COD concentrations below the imposed limit of 50 g m-3 whilst the
greatest predictions of effluent COD concentration were as high as 900 g m-3
and concentrated within 25% of the search space. This may explain the high GPR
errors for these variables due to these high magnitude outliers resulting in
high MAPE errors relative to most of the training data.

The maximum limit for the effluent TN concentration (10 g m-3) results in
infeasible process designs at low UASB volumes for high influent COD
concentrations as well as at high UASB volumes for lower influent COD
concentrations. The greatest effluent TN concentrations were observed for the
UCT process with low UASB volumes and high influent COD concentrations.

The TP effluent quality limit of 5 g m-3 enforced infeasible regions again at
low UASB volumes and high influent COD concentrations as well as at high UASB
volumes, particularly for low COD concentrations across all 4 BNR pathways,
yet impacting feasible designs at high UASB volumes and high COD
concentrations for A2O, Bardenpho, and Johannesburg configurations. The
greatest violations of the TP effluent quality limit were observable at high
UASB volumes and low influent COD concentrations for the 3 aforementioned BNR
process configurations.

The classification NN results are shown for the 4 BNR pathways in Figure 13.
The NN demonstrates a good separation between feasible and infeasible designs
at logit predictions of 0 as shown by the separating hyperplane. Since all of
the samples converged for the UCT pathway, the logit predictions are greater
than 0 across the entire input space. The classification NNs were validated
based on precision and recall scores which were both 1 across Bardenpho,
Johannesburg, and UCT pathways, whilst A2O exhibited a precision of 0.88 (due
to the incorrectly classified non-converged sample within the testing set
visible in Figure 13) and a recall of 1.

Figure 13: Surrogate modelling of the feasible region (shown by the dashed
lines) constrained by an effluent COD limit of 50 g m-1, an effluent TN limit
of 10 g m-1, and effluent TP limit of 5 g m-1, and a probability of simulation
convergence greater than 0.5 (corresponding to a logit prediction of 0).
Effluent COD, TN, and TP were modelled using GPR whilst the logit predictions
were from classification NNs. Note that the UCT simulations converged for all
samples therefore there is no region of infeasibility due to simulation
convergence failures. Converged samples are shown as black whilst non-
converged samples are shown as red. Additionally, training samples are shown
as filled whilst testing samples are shown as unfilled.

Figure 14 shows the GPR response surfaces for the objective functions
enumerating the nutrient quality of the recovered digestate, the biogas
recovery, and the aeration requirement of different process designs within the
superstructure. Unconstrained by the feasibility constraints, the maximum
nutrient quality was achieved at low UASB volumes and low influent COD
concentrations for the Bardenpho process. Maximum biogas recovery was realised
at high UASB volumes and high influent COD concentrations whilst the minimum
aeration flow was required for the upper half of UASB volumes over a range of
influent COD concentrations expect the greatest values.

Figure 14: Objective function GPR response surfaces.

Figure 15 shows the standard deviations in the GPR predictions from Figure 14.
The interpolation property of GPR models is observable by the negligible
standard deviations close to the training samples with increasing
uncertainties with increasing distance from samples. Note that the non-
converged samples were not used to train the GPR models such that the
uncertainty in the GPR predictions was unaffected by proximity to these
samples. Because the same training data was used for each GPR model, the shape
of the standard deviation response surfaces are all similar, only varying in
the amplitude of uncertainties based on the magnitude of the corresponding
variable predictions.

Figure 15: GPR standard deviation response surfaces. The training data are
superimposed where non-converged training samples (used for classification
model training but not used to train the regression models) are shown in red.

###  4.3 Multi-objective optimisation results

Figure 16 shows the results of the superstructure optimisation to maximise the
combined N and P nutrient quality of the recovered digestate at 61.2 $\pm$ 2.5
kg t-1 (where the uncertainty was obtained from the standard deviation in the
GPR predictions based on the 95% confidence interval). The optimal flowsheet
for nutrient recovery implemented the Bardenpho BNR process with a 143 m3 UASB
and was realised at an influent COD of 2,000 g m-3 as shown by the red star
(Figure 16). The amount of recovered biogas with this flowsheet is 165 m3d-1
whilst the aeration requirement is 41,340 m3d-1. It can be seen from Figure 16
that none of the constraints are active since the flowsheet produces an
effluent with concentrations of COD, TN, and TP of 20.6 g m-3, 3.4 g m-3, and
1.6 g m-3, respectively.

Figure 16 also shows the results of the superstructure optimisation to
maximise the total biogas recovery as 1,500 $\pm$ 103 m3 d-1. The maximal
biogas recovery was obtained using the Johannesburg BNR process with a 388 m3
UASB and was realised at an influent COD of 5,616 g m-3 as shown by the red
star (Figure 16). The nutrient quality of the recovered digestate using this
flowsheet is 26.5 kg t-1 whilst the aeration requirement is 18,920m3d-1. None
of the effluent constraints were active for this solution since the flowsheet
produced an effluent with concentrations of TN and TP of 4.9 g m-3, and 0.2 g
m-3, respectively, whilst all of the COD was predicted to have been recovered.

Figure 16 also shows the results of the superstructure optimisation to
minimise the total aeration requirement as 15,790 $\pm$ 8,080 m3 d-1. This
solution was obtained using the UCT BNR process with a 255 m3 UASB and was
realised at an influent COD of 2,598 g m-3 as shown by the red star (Figure
16). The nutrient quality of the recovered digestate and the amount of
recovered biogas using this flowsheet was 39.7 kg t-1 and 584 m3d-1,
respectively. The effluent constraints for this solution had concentrations of
COD, TN, and TP of 12.8 g m-3, 9.8 g m-3, and only 5.0 g m-3, respectively.
The TP effluent constraint was therefore active at this solution which can
been seen in Figure 16.

A Pareto optimal flowsheet was determined to recover a digestate with an N and
P nutrient quality of 35.9 $\pm$ 6.5 kg t-1 alongside 1,070 m3d-1 of biogas,
with an aeration requirement of 31,570 m3d-1. This Pareto optimal flowsheet
implemented to UCT process with a UASB volume of 360 m3d-1, and operated
optimally for an influent COD concentration of 4,871 g m-3 (Figure 16). The
effluent from the Pareto optimal flowsheet had COD, TN, and TP concentrations
of 0 g m-3, 2.4 g m-3, and 0 g m-3, respectively.

Figure 16: Optimisation solutions for 4 decision criteria: maximised nutrient
quality (upper left); maximised biogas recovery (upper right); minimised
aeration cost (lower left); and a Pareto optimal solution to maximise nutrient
quality whilst ensuring a biogas recovery of at least 50 % the optimal (lower
right). The optimised UASB volumes and corresponding influent COD
concentrations are shown by the red starts whilst the optimal BNR pathway
appears in the title above each response surface. The agglomeration of the
effluent quality feasible regions and the simulation convergence-based
feasible region is constrained within the black shaded regions.

It is also worth noting the different influent COD concentrations realised for
the optimal solutions in Figure 16. Specifically, the optimal nutrient
recovery was realised at the lower bound of 2,000 g m-3 COD whilst the maximum
biogas recovery was obtained for a COD concentration nearer the upper bound at
5,616 g m-3. Therefore, nutrient recovery could be the focus of the process
system at lower influent COD concentrations, whilst biogas recovery is
favoured at higher influent COD concentrations. To this end, a flexible
resource recovery process, in which Bardenpho and Johannesburg BNR
configurations were configured, could enable the operation of the resource
recovery process to vary temporally based on the influent COD. However, this
would require optimisation of the process system operations in response to
real-time sensor data and future projections, for example harnessing model
predictive control.

Figure 17 shows the trade-offs between the competing objective criteria where
each single objective solution is visualised as well as the Pareto optimal
solution. The maximised nutrient quality of the recovered digestate was 61.2
kg t-1 whilst the most sub-optimal nutrient quality existed for the solution
to maximise biogas recovery at 26.5 kg t-1 (a reduction of over 55 %).
Similarly, the maximum biogas recovery was 1,500 m3 d-1 whilst the lowest
biogas recovery existed for the flowsheet to maximise the nutrient recovery at
164.5 m3d-1 or just 11 % of the maximum biogas recovery potential. Finally,
the minimised aeration requirement was 15,790 m3 d-1 whilst the greatest
aeration requirement was almost 3 times as high and existed for the maximal
nutrient recovery solution at 41,340 m3d-1. The Pareto optimal solution
achieved a trade-off between the 3 objectives with 59 % of the maximised
nutrient recovery, 71 % of the maximum biogas recovery, necessitating 76 % of
the aeration requirement to recover the maximum nutrient quality digestate.

Figure 17: Multi-objective optimisation trade-offs. Values are normalised to
the radial axes where solutions furthest from the origin are optimal and
solutions nearest to the origin are the most sub-optimal trade-offs.

###  4.4 Optimisation under uncertainty results

Solutions to the optimisation under influent COD uncertainties determined the
optimum BNR configuration and UASB volumes. The solution to the 3 objectives
to maximise nutrient quality, biogas recovery, and minimise aeration cost
converged on the same solution implementing a 274 m3 UASB followed by the UCT
BNR process. Figure 18 thereby shows the GPR predictions for the optimal
process which recovers a digestate with an expected nutrient quality of 38.6
kg t-1 (63 % of the maximised nutrient quality without any uncertainty
considerations), 580 m3 d-1 expected biogas recovery (39 % of the maximum
value found previously), and necessitating an expected aeration flow of 40,360
m3 d-1 which was over 150 % greater than the minimised value without any
uncertainty considerations. Specifically, Figure 18 shows the responses of the
optimal process design to the different realisations of influent COD between
2,000 g m-3 and 6,000 g m-3.

The expected effluent concentrations were constrained to be less than the
effluent limits (50 g m-3, 10 g m-3, and 10 g m-3 for COD, TN, and TP,
respectively). Figure 18 shows that the effluent constraint for the COD
concentration was below its specified limit with an expected value of 48.1 g
m-3. Additionally, it can be seen that the effluent COD is expected to
transgress the limits at influent COD concentrations greater than 5,000 g m-3.

Similarly, despite the expected effluent TN concentration of 7.6 g m-3 being
below the effluent limit of 10 g m-3, the GPR predictive mean violates this
limit at influent COD concentrations above 5,000 g m-3. Similarly the expected
effluent TP concentration was 2.4 g m-3 whilst the predictive mean function
over the realisations, despite not violating the effluent limit of 10 g m-3
over all realisations of the influent COD, the 95 % upper confidence bound
violates this limit for extreme influent COD concentrations. In fact,
infringement of the effluent limits by the 95 % upper confidence bound is
common across COD, TN, and TP concentrations.

Figure 18 also shows the predictive GPR model for the recovered digestate
nutrient quality and the uncertainty in these predictions dependent on
realisations of the influent COD concentration. Also shown is the expected
value across all realisations of the influent COD, accounting for the
probability of these concentrations being realised. Generally, the nutrient
quality appears to greatest for lower influent COD concentrations, whilst the
nutrient quality at greater influent COD concentrations is reduced. In
particular, the uncertainty in the nutrient quality peaks at extreme influent
COD concentrations of 2,000 g m-3 and 6,000 g m-3. In fact, this
characteristic of high uncertainty at the bounds of the influent COD
realisations is common across all of the GPR models due to the dependence on
the distribution of training data sampled from the process simulation
software. The GPR modelling uncertainties could thereby be reduced by
increasing the number of static samples generated in initial sampling
strategies or by adopting an adaptive sampling approach wherein additional
samples are iteratively evaluated from the process simulator so as to increase
the modelling fidelity at sparsely sampled locations.

Figure 18: Operational profiles of the solution to the optimisation problem
accounting for uncertainty in the influent chemical oxygen demand (COD). Also
shown are the predictive modelling uncertainties (as the 95 % confidence
interval) from Gaussian process regression models. TN: total nitrogen, TP:
total phosphorous.

##  5 Conclusions

An object-orientated methodology was developed for surrogate modelling and
DFO. GP and NN surrogate models were developed for both regression and
classification applications. Explicit mathematical formulations of these
surrogate models were validated against existing machine learning
implementations and embedded within generalised optimisation formulations.
Specifically, GP models were formulated as NLPs whilst NN models were
formulated as MILPs or NLPs depending on whether the activation functions were
piecewise linear or nonlinear, respectively. These abstracted optimisation
formulations enabled plug-and-play surrogate-based optimisation modelling. A
primary contribution of this work was the development of classification
surrogate models to enable feasibility constraint formulations addressing
uncertainties pertaining to simulation convergence failures.

Adaptive sampling algorithms based on GPs and Delaunay triangulation were
developed to improve the exploration and exploitation of the search space. GP-
based adaptive sampling methods enabled exploration via maximisation of the GP
uncertainty or exploitation via the modified EI, both formulated as NLPs.
Delaunay triangulation methods were formulated as MILPs to enable heuristic-
based exploration and exploitation of the search space. The plug-and-play
feature of the developed optimisation formulations enabled online feasibility
constraints to be plugged-in to adaptive sampling formulations and increase
the efficiency of sampling in the feasible region. Finally, the heuristic-
based adaptive sampling MILP was adapted within a data-based direct-search
optimisation algorithm to optimise black box problems. The adoption of the
object-orientated paradigm means that these methods can be harnessed in a
plug-and-play approach to increase the global optimisation and sampling
efficiencies of a wide range of simulation-based optimisation applications.

This work demonstrated the applicability of the methodology developed to the
optimisation of process systems to recover resources from food and beverage
processing wastewaters (FPWW). Specifically, a superstructure optimisation
methodology was developed to optimise process systems to recover carbon,
nitrogen, and phosphorous resources from a typical brewery wastewater. A
superstructure embedding a high-rate anaerobic digestion reaction upstream of
4 different BNR configurations (A2O, Bardenpho, Johannesburg, and UCT) was
implemented in GPS-X process simulation software so as to harness high-
fidelity data within the decision making framework. This work therefore
contributed a practical simulation-based process synthesis application of a
black box mixed-integer nonlinear programming framework embedding both GP and
NN surrogate model formulations. To achieve this, GPR surrogate models were
harnessed to represent resource recovery pathways with interpretable
uncertainty quantification. Additionally, classification NN surrogate models
were used to model the feasible region whilst addressing uncertainties from
non-converged black box simulation data.

Multi-objective optimisation addressed the trade-offs between 3 competing
decision criteria: one to maximise the total recovered biogas; one to maximise
the nutrient quality of a recoverable digestate; and one to minimise the
aeration requirement as a proxy to the total operational cost. As such, this
work addressed the challenge of holistic process synthesis which currently
constrains the deployment of resource recovery process to wastewater systems.
Optimal nutrient recovery was realised at low COD compositions, whilst maximum
biogas recovery was achieved at COD compositions near to the upper bound cited
for typical brewery effluents.

Uncertainties in the brewery wastewater COD composition were accounted for by
implementing a stochastic programming reformulation of the optimisation
problem, thereby enabling the process synthesis to account for the inherent
uncertainties in wastewater characterisation. This work therefore demonstrates
the applicability of optimisation under uncertainty methods to resource
recovery from FPWW thereby contributing to the literature which focuses
primarily on applications to organic solid wastes. The expected effluent
compositions were enforced to be below maximum concentration limits
accordingly. The optimisation also enforced robustness to process feasibility
constraints derived from convergence of the process simulator. The solutions
to the stochastic programming formulations maximised the expected nutrient
recovery or biogas recovery, or minimised the expected aeration cost. The same
solution, a 274 m3 UASB upstream of a UCT BNR configuration, was determined
for all 3 objectives. This was suspected due to the difficulty of satisfying
the expected effluent constraints over realisations of the influent COD
uncertainties. The modelling framework enables machine learning models of the
optimal process design to be harnessed for further operations optimisation.
Specifically, GPR models provided predictions of process performance in
response to variations in the influent COD composition in addition to the
uncertainties in these predictions.

###  5.1 Recommendations for future work

Flexible process optimisation was recommended to determine flexible resource
recovery processes in which the BNR pathway can be adjusted in response to
temporal variations in the wastewater composition. Another recommendation was
to integrate the uncertainties in the GPR surrogate models within the
optimisation-based decision-making process. One strategy to achieve this is
via minimising the upper confidence bound which can be calculated by adding
factors of the GPR standard deviations to the GPR predictive mean. For
example, the 95 % upper confidence bound can be obtained by adding 1.96
standard deviations to the mean predictions provided by the GPR modelling
objects. This would ensure that obtained solutions simultaneously minimised
the objective function and the uncertainty in solutions. The upper confidence
bound could also be used within the constraint formulations to ensure the the
optimisation solutions were robust to the GPR modelling uncertainties that
arise due to interpolation between the training data.

There exist multiple research directions to improve the capabilities of the
developed methodology. These include future work to increase the number of
surrogate models available within the object-orientated framework. One example
to achieve this is to consider enabling the selection of more kernel functions
within GPC models, particularly a linear kernel which would reduce the
complexity of GPC-based NLP formulations. Another example would be to enable
recurrent NN structures to improve nonlinear deep learning capabilities. A
third example is to incorporate other surrogate models into the framework for
regression and classification: support vector machines, random forests, and
decision trees have been highlighted as additional alternatives.

Future work could also focus on improving the implementations of the GP and NN
surrogate models. This could be achieved by enabling optimisation of the order
of the GPR polynomial kernel. Another option is to enable the optimisation of
NN parameters such as number of layers, nodes in each layer, and activation
functions [112]. This could further be adapted to implement algorithms to
optimise the selection of different surrogate models and their respective
model structure as in [63]. The surrogate modelling capabilities could be
further improved by implementing further regularisation methods such as L1 and
dropout regularisation for NNs as well as early termination prior to
overfitting. Cross-validation methods are another method to validate and
quantify overfitting of surrogate models. Additionally, analysis of the effect
of multiple NN outputs compared to multiple single-output NNs on the resulting
NN-based optimisation formulations should be explored. Other adaptive sampling
formulations and acquisition functions could also improve the existing
codebase’s capabilities. [43] provide a good discussion on the variants of GP
models using different correlation and mean functions during specification of
the GP prior. Additionally, [43] provide a good discussion on the
computational aspects on GP modelling pertaining to handling higher
dimensional problems and the non-convexity of the maximum likelihood
estimator.

Finally, future work should implement the open-source methodology on more
applications. These could be programmatic testing applications to inform
further programmable improvements of the codebase and its wider usability. In
this regard, increased documentation and demonstrative applications would be
beneficial. Application of the methodology to real world engineering problems
is also particularly important to realise the impact of the developed
methodology in practical operating environments. To this end, it would be
beneficial to combine these data-driven surrogate models with well-established
mechanistic models and expert knowledge. Such hybrid models would improve
transparency, interpretability and extrapolation capabilities thereby gaining
the trust of practical operational engineers with a lack of computer science
skills. Finally, such hybrid models should leverage the increasing volumes of
available sensor data. This would remove a layer of abstraction incurred from
the reliance on simulation models and thereby further increase the
transparency of black box modelling frameworks.

## References

  * [1] S. Joe Qin and Leo H. Chiang.  Advances and opportunities in machine learning for process data analytics.  Computers & Chemical Engineering, 126:465–473, July 2019.
  * [2] E N Pistikopoulos, Ana Barbosa-Povoa, Jay H Lee, Ruth Misener, Alexander Mitsos, G V Reklaitis, V Venkatasubramanian, Fengqi You, and Rafiqul Gani.  Process systems engineering – The generation next?  Computers & Chemical Engineering, 147:107252, April 2021.
  * [3] Artur M. Schweidtmann, Erik Esche, Asja Fischer, Marius Kloft, Jens-Uwe Repke, Sebastian Sager, and Alexander Mitsos.  Machine Learning in Chemical Engineering: A Perspective.  Chemie Ingenieur Technik, 93(12):2029–2039, 2021.
  * [4] Mariane Yvonne Schneider, Ward Quaghebeur, Sina Borzooei, Andreas Froemelt, Feiyi Li, Ramesh Saagi, Matthew J. Wade, Jun-Jie Zhu, and Elena Torfs.  Hybrid modelling of water resource recovery facilities: status and opportunities.  Water Science and Technology, 85(9):2503–2524, April 2022.
  * [5] Francesco Cherubini, Gerfried Jungmeier, Maria Wellisch, Thomas Willke, Ioannis Skiadas, René Van Ree, and Ed de Jong.  Toward a common classification approach for biorefinery systems.  Biofuels, Bioproducts and Biorefining, 3(5):534–546, 2009.
  * [6] Jian Gong and Fengqi You.  Sustainable design and synthesis of energy systems.  Current Opinion in Chemical Engineering, 10:77–86, November 2015\.
  * [7] Jian Gong and Fengqi You.  Value-Added Chemicals from Microalgae: Greener, More Economical, or Both?  ACS Sustainable Chemistry & Engineering, 3(1):82–96, January 2015\.
  * [8] Atul Kumar and S. R. Samadder.  A review on technological options of waste to energy for effective management of municipal solid waste.  Waste Management, 69:407–422, November 2017.
  * [9] Robert Pott, Madelyn Johnstone-Robertson, Bernelle Verster, Shilpa Rumjeet, Lefa Nkadimeng, Tayana Raper, Sharon Rademeyer, and Susan T. L. Harrison.  Wastewater Biorefineries: Integrating Water Treatment and Value Recovery.  In Walter Leal Filho and Dinesh Surroop, editors, The Nexus: Energy, Environment and Climate Change, Green Energy and Technology, pages 289–302. Springer International Publishing, Cham, 2018.
  * [10] Yi-di Chen, Shih-Hsin Ho, Dillirani Nagarajan, Nan-qi Ren, and Jo-Shu Chang.  Waste biorefineries — integrating anaerobic digestion and microalgae cultivation for bioenergy production.  Current Opinion in Biotechnology, 50:101–110, April 2018.
  * [11] J. Rajesh Banu, S. Kavitha, R. Yukesh Kannah, Rahul R. Bhosale, and Gopalakrishnan Kumar.  Industrial wastewater to biohydrogen: Possibilities towards successful biorefinery route.  Bioresource Technology, 298:122378, February 2020.
  * [12] Luca Alibardi, Thomas F. Astrup, Fabiano Asunis, William P. Clarke, Giorgia De Gioannis, Paolo Dessì, Piet N. L. Lens, Maria Cristina Lavagnolo, Lidia Lombardi, Aldo Muntoni, Alberto Pivato, Alessandra Polettini, Raffaella Pomi, Andreina Rossi, Alessandro Spagni, and Daniela Spiga.  Organic waste biorefineries: Looking towards implementation.  Waste Management, 114:274–286, August 2020.
  * [13] Hui Yi Leong, Chih-Kai Chang, Kuan Shiong Khoo, Kit Wayne Chew, Shir Reen Chia, Jun Wei Lim, Jo-Shu Chang, and Pau Loke Show.  Waste biorefinery towards a sustainable circular bioeconomy: a solution to global issues.  Biotechnology for Biofuels, 14(1):87, April 2021.
  * [14] C. Puchongkawarin, C. Gomez-Mont, D. C. Stuckey, and B. Chachuat.  Optimization-based methodology for the development of wastewater facilities for energy and nutrient recovery.  Chemosphere, 140:150–158, December 2015.
  * [15] A. Alva-Argáez, A. C. Kokossis, and R. Smith.  Wastewater minimisation of industrial systems using an integrated approach.  Computers & Chemical Engineering, 22:S741–S744, March 1998.
  * [16] Hector Yeomans and Ignacio E. Grossmann.  A systematic modeling framework of superstructure optimization in process synthesis.  Computers & Chemical Engineering, 23(6):709–731, June 1999.
  * [17] Qi Chen and I. E. Grossmann.  Recent Developments and Challenges in Optimization-Based Process Synthesis.  Annual Review of Chemical and Biomolecular Engineering, 8:249–283, June 2017.
  * [18] Elvis Ahmetović and Ignacio E. Grossmann.  Global superstructure optimization for the design of integrated process water networks.  AIChE Journal, 57(2):434–457, 2011.
  * [19] Hande Bozkurt, Alberto Quaglia, Krist V. Gernaey, and Gürkan Sin.  A mathematical programming framework for early stage design of wastewater treatment plants.  Environmental Modelling & Software, 64:164–176, February 2015\.
  * [20] Hande Bozkurt, Mark C. M. van Loosdrecht, Krist V. Gernaey, and Gürkan Sin.  Optimal WWTP process selection for treatment of domestic wastewater – A realistic full-scale retrofitting study.  Chemical Engineering Journal, 286:447–458, February 2016.
  * [21] A. Castillo, P. Cheali, V. Gómez, J. Comas, M. Poch, and G. Sin.  An integrated knowledge-based and optimization tool for the sustainable selection of wastewater treatment process concepts.  Environmental Modelling & Software, 84:177–192, October 2016.
  * [22] Juan I. Padrón-Páez, Sofía De-León Almaraz, and Alicia Román-Martínez.  Sustainable wastewater treatment plants design through multiobjective optimization.  Computers & Chemical Engineering, 140:106850, September 2020.
  * [23] Resul Al, Chitta Ranjan Behera, Krist V. Gernaey, and Gürkan Sin.  Stochastic simulation-based superstructure optimization framework for process synthesis and design under uncertainty.  Computers & Chemical Engineering, 143:107118, December 2020.
  * [24] Živko Južnič-Zonta, Albert Guisasola, and Juan Antonio Baeza.  Smart-Plant Decision Support System (SP-DSS): Defining a multi-criteria decision-making framework for the selection of WWTP configurations with resource recovery.  Journal of Cleaner Production, 367:132873, September 2022.
  * [25] Alex Durkin, Tim Finnigan, Rob Johnson, John Kazer, Jianfeng Yu, David Stuckey, and Miao Guo.  Can closed-loop microbial protein provide sustainable protein security against the hunger pandemic?  Current Research in Biotechnology, September 2022.
  * [26] D. J. Batstone, J. Keller, I. Angelidaki, S. V. Kalyuzhnyi, S. G. Pavlostathis, A. Rozzi, W. T. M. Sanders, H. Siegrist, and V. A. Vavilin.  The IWA Anaerobic Digestion Model No 1 (ADM1).  Water Science and Technology, 45(10):65–73, May 2002.
  * [27] M. Henze, W. Gujer, T. Mino, and M. van Loosedrecht.  Activated Sludge Models ASM1, ASM2, ASM2d and ASM3.  IWA Publishing, October 2006.
  * [28] Kimberly Solon, Eveline I. P. Volcke, Mathieu Spérandio, and Mark C. M. van Loosdrecht.  Resource recovery and wastewater treatment modelling.  Environmental Science: Water Research & Technology, 5(4):631–642, March 2019.
  * [29] Saman Razavi, Bryan A. Tolson, and Donald H. Burn.  Review of surrogate modeling in water resources.  Water Resources Research, 48(7), 2012.
  * [30] Jerome Sacks, William J. Welch, Toby J. Mitchell, and Henry P. Wynn.  Design and Analysis of Computer Experiments.  Statistical Science, 4(4):409–423, November 1989.
  * [31] Satyajith Amaran, Nikolaos V. Sahinidis, Bikram Sharda, and Scott J. Bury.  Simulation optimization: a review of algorithms and applications.  Annals of Operations Research, 240(1):351–380, May 2016.
  * [32] Antonia Hadjimichael, Joaquim Comas, and Lluís Corominas.  Do machine learning methods used in data mining enhance the potential of decision support systems? A review for the urban water sector.  AI Communications, 29(6):747–756, January 2016.
  * [33] Ll. Corominas, M. Garrido-Baserba, K. Villez, G. Olsson, U. Cortés, and M. Poch.  Transforming data into knowledge for improved wastewater treatment operation: A critical review of techniques.  Environmental Modelling & Software, 106:89–103, August 2018.
  * [34] Kathryn B. Newhart, Ryan W. Holloway, Amanda S. Hering, and Tzahi Y. Cath.  Data-driven performance analyses of wastewater treatment plants: A review.  Water Research, 157:498–513, June 2019.
  * [35] Mengyuan Zhu, Jiawei Wang, Xiao Yang, Yu Zhang, Linyu Zhang, Hongqiang Ren, Bing Wu, and Lin Ye.  A review of the application of machine learning in water quality evaluation.  Eco-Environment & Health, 1(2):107–116, June 2022.
  * [36] Maarten R. Dobbelaere, Pieter P. Plehiers, Ruben Van de Vijver, Christian V. Stevens, and Kevin M. Van Geem.  Machine Learning in Chemical Engineering: Strengths, Weaknesses, Opportunities, and Threats.  Engineering, 7(9):1201–1211, September 2021.
  * [37] William Bradley, Jinhyeun Kim, Zachary Kilwein, Logan Blakely, Michael Eydenberg, Jordan Jalvin, Carl Laird, and Fani Boukouvala.  Perspectives on the integration between first-principles and data-driven modeling.  Computers & Chemical Engineering, 166:107898, October 2022.
  * [38] David M. Himmelblau.  Accounts of Experiences in the Application of Artificial Neural Networks in Chemical Engineering.  Industrial & Engineering Chemistry Research, 47(16):5782–5796, August 2008.
  * [39] Xiupeng Wei.  Modeling and optimization of wastewater treatment process with a data-driven approach.  PhD thesis, University of Iowa, 2013.
  * [40] Andrew R. Conn, Katya Scheinberg, and Luis N. Vicente.  Introduction to derivative-free optimization.  SIAM, 2009.
  * [41] Charles Audet and Warren Hare.  Derivative-Free and Blackbox Optimization.  Springer Series in Operations Research and Financial Engineering. Springer International Publishing, Cham, 2017.
  * [42] Luis Miguel Rios and Nikolaos V. Sahinidis.  Derivative-free optimization: a review of algorithms and comparison of software implementations.  Journal of Global Optimization, 56(3):1247–1293, July 2013.
  * [43] Atharv Bhosekar and Marianthi Ierapetritou.  Advances in surrogate based modeling, feasibility analysis, and optimization: A review.  Computers & Chemical Engineering, 108:250–267, January 2018.
  * [44] Stéphane Alarie, Charles Audet, Aïmen E. Gheribi, Michael Kokkolaras, and Sébastien Le Digabel.  Two decades of blackbox optimization applications.  EURO Journal on Computational Optimization, 9:100011, January 2021\.
  * [45] Peter I. Frazier.  A Tutorial on Bayesian Optimization, July 2018.
  * [46] Sun Hye Kim and Fani Boukouvala.  Surrogate-based optimization for mixed-integer nonlinear problems.  Computers & Chemical Engineering, 140:106847, September 2020.
  * [47] Fani Boukouvala, Ruth Misener, and Christodoulos A. Floudas.  Global optimization advances in Mixed-Integer Nonlinear Programming, MINLP, and Constrained Derivative-Free Optimization, CDFO.  European Journal of Operational Research, 252(3):701–727, August 2016.
  * [48] Fani Boukouvala, M. M. Faruque Hasan, and Christodoulos A. Floudas.  Global optimization of general constrained grey-box models: new method and its application to constrained PDEs for pressure swing adsorption.  Journal of Global Optimization, 67(1):3–42, January 2017.
  * [49] Artur M. Schweidtmann and Alexander Mitsos.  Deterministic Global Optimization with Artificial Neural Networks Embedded.  Journal of Optimization Theory and Applications, 180(3):925–948, March 2019.
  * [50] Donato Maragno, Holly Wiberg, Dimitris Bertsimas, S. Ilker Birbil, Dick den Hertog, and Adejuyigbe Fajemisin.  Mixed-Integer Optimization with Constraint Learning, November 2021\.
  * [51] Francesco Ceccon, Jordan Jalving, Joshua Haddad, Alexander Thebelt, Calvin Tsay, Carl D. Laird, and Ruth Misener.  OMLT: Optimization & Machine Learning Toolkit, February 2022\.  arXiv:2202.02414 [cs, math, stat].
  * [52] Michael L. Bynum, Gabriel A. Hackebeil, William E. Hart, Carl D. Laird, Bethany L. Nicholson, John D. Siirola, Jean-Paul Watson, and David L. Woodruff.  Pyomo — Optimization Modeling in Python, volume 67.  Springer Nature, 2021.
  * [53] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [54] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  PyTorch: An Imperative Style, High-Performance Deep Learning Library.  In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [55] Bjarne Grimstad and Henrik Andersson.  ReLU Networks as Surrogate Models in Mixed-Integer Linear Programs.  Computers & Chemical Engineering, 131:106580, December 2019.
  * [56] Vincent Tjeng, Kai Xiao, and Russ Tedrake.  Evaluating Robustness of Neural Networks with Mixed Integer Programming, February 2019.
  * [57] José A. Caballero and Ignacio E. Grossmann.  An algorithm for the use of surrogate models in modular flowsheet optimization.  AIChE Journal, 54(10):2633–2650, 2008.
  * [58] Carlos A. Henao and Christos T. Maravelias.  Surrogate-based superstructure optimization framework.  AIChE Journal, 57(5):1216–1232, 2011.
  * [59] Kenneth Holmström and Marcus M. Edvall.  The TOMLAB Optimization Environment.  In Josef Kallrath, editor, Modeling Languages in Mathematical Optimization, Applied Optimization, pages 369–376. Springer US, Boston, MA, 2004.
  * [60] Philip E. Gill, Walter Murray, and Michael A. Saunders.  SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization.  SIAM Journal on Optimization, 12(4):979–1006, January 2002.
  * [61] Michael R. Bussieck and Alex Meeraus.  General Algebraic Modeling System (GAMS).  In Josef Kallrath, editor, Modeling Languages in Mathematical Optimization, Applied Optimization, pages 137–157. Springer US, Boston, MA, 2004.
  * [62] Ignacio E Grossmann, Jagadisan Viswanathan, Aldo Vecchietti, Ramesh Raman, and Erwin Kalvelagen.  GAMS/DICOPT: A Discrete Continuous Optimization Package, 2002.
  * [63] Fani Boukouvala and Christodoulos A. Floudas.  ARGONAUT: AlgoRithms for Global Optimization of coNstrAined grey-box compUTational problems.  Optimization Letters, 11(5):895–913, June 2017.
  * [64] Burcu Beykal, Melis Onel, Onur Onel, and Efstratios N. Pistikopoulos.  A data-driven optimization algorithm for differential algebraic equations with numerical infeasibilities.  AIChE Journal, 66(10):e16657, 2020.
  * [65] Ruth Misener and Christodoulos A. Floudas.  ANTIGONE: Algorithms for coNTinuous / Integer Global Optimization of Nonlinear Equations.  Journal of Global Optimization, 59(2):503–526, July 2014.
  * [66] Dirk Gorissen, Ivo Couckuyt, Piet Demeester, Tom Dhaene, and Karel Crombecq.  A Surrogate Modeling and Adaptive Sampling Toolbox for Computer Based Design.  The Journal of Machine Learning Research, 11:2051–2055, August 2010\.
  * [67] Mohamed Amine Bouhlel, John T. Hwang, Nathalie Bartoli, Rémi Lafage, Joseph Morlier, and Joaquim R. R. A. Martins.  A Python surrogate modeling framework with derivatives.  Advances in Engineering Software, 135:102662, September 2019.
  * [68] Alison Cozad, Nikolaos V. Sahinidis, and David C. Miller.  Learning surrogate models for simulation-based optimization.  AIChE Journal, 60(6):2211–2227, 2014.
  * [69] Charles Audet, Sébastien Le Digabel, Viviane Rochon Montplaisir, and Christophe Tribes.  NOMAD version 4: Nonlinear optimization with the MADS algorithm, May 2021.
  * [70] Fani Boukouvala and Marianthi G. Ierapetritou.  Surrogate-Based Optimization of Expensive Flowsheet Modeling for Continuous Pharmaceutical Manufacturing.  Journal of Pharmaceutical Innovation, 8(2):131–145, June 2013.
  * [71] J. Edward Jackson.  A User’s Guide to Principal Components.  Wiley, March 1991.
  * [72] Matthias Schonlau and William J. Welch.  Screening the Input Variables to a Computer Model Via Analysis of Variance and Visualization.  In Angela Dean and Susan Lewis, editors, Screening: Methods for Experimentation in Industry, Drug Discovery, and Genetics, pages 308–327. Springer, New York, NY, 2006.
  * [73] I. M Sobol.  Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates.  Mathematics and Computers in Simulation, 55(1):271–280, February 2001.
  * [74] Alexander I. J. Forrester and Andy J. Keane.  Recent advances in surrogate-based optimization.  Progress in Aerospace Sciences, 45(1):50–79, January 2009.
  * [75] M. D. McKay, R. J. Beckman, and W. J. Conover.  A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code.  Technometrics, 21(2), 1979.
  * [76] Jayant R. Kalagnanam and Urmila M. Diwekar.  An Efficient Sampling Technique for Off-line Quality Control.  Technometrics, 39(3), August 1997.
  * [77] I. M Sobol’.  On the distribution of points in a cube and the approximate evaluation of integrals.  USSR Computational Mathematics and Mathematical Physics, 7(4):86–112, January 1967.
  * [78] Sushant S. Garud, Iftekhar A. Karimi, and Markus Kraft.  Design of computer experiments: A review.  Computers & Chemical Engineering, 106:71–95, November 2017.
  * [79] M. E. Johnson, L. M. Moore, and D. Ylvisaker.  Minimax and maximin distance designs.  Journal of Statistical Planning and Inference, 26(2):131–148, October 1990.
  * [80] Donald R. Jones, Matthias Schonlau, and William J. Welch.  Efficient Global Optimization of Expensive Black-Box Functions.  Journal of Global Optimization, 13(4):455–492, December 1998.
  * [81] Fani Boukouvala and Marianthi G. Ierapetritou.  Derivative-free optimization for expensive constrained problems using a novel expected improvement objective function.  AIChE Journal, 60(7):2462–2474, 2014.
  * [82] Sushant Suhas Garud, I. A. Karimi, and Markus Kraft.  Smart Sampling Algorithm for Surrogate Model Development.  Computers & Chemical Engineering, 96:103–114, January 2017.
  * [83] B. Delaunay.  Sur la sphère vide.  Bulletin de l’Académie des Sciences de l’URSS. Classe des sciences mathématiques et na, 1934(6):793–800, 1934.
  * [84] Ping Jiang, Yahui Zhang, Qi Zhou, Xinyu Shao, Jiexiang Hu, and Leshi Shu.  An adaptive sampling strategy for Kriging metamodel based on Delaunay triangulation and TOPSIS.  Applied Intelligence, 48(6):1644–1656, June 2018.
  * [85] Sushant S. Garud, Nivethitha Mariappan, and Iftekhar A. Karimi.  Surrogate-based black-box optimisation via domain exploration and smart placement.  Computers & Chemical Engineering, 130:106567, November 2019.
  * [86] Pooriya Beyhaghi, Daniele Cavaglieri, and Thomas Bewley.  Delaunay-based derivative-free optimization via global surrogates, part I: linear constraints.  Journal of Global Optimization, 66(3):331–382, November 2016.
  * [87] Eddie Davis and Marianthi Ierapetritou.  A centroid-based sampling strategy for kriging global modeling and optimization.  AIChE Journal, 56(1):220–240, 2010.
  * [88] Y. Wu, L. Ozdamar, and A. Kumar.  TRIOPT: a triangulation-based partitioning algorithm for global optimization.  Journal of Computational and Applied Mathematics, 177(1):35–53, May 2005.
  * [89] Bernd Gärtner and Michael Hoﬀmann.  Computational Geometry Lecture Notes HS 2013, 2013.
  * [90] Andres Gonzalez-Garay and Gonzalo Guillen-Gosalbez.  SUSCAPE: A framework for the optimal design of SUStainable ChemicAl ProcEsses incorporating data envelopment analysis.  Chemical Engineering Research and Design, 137:246–264, September 2018.
  * [91] Raul Yondo, Esther Andrés, and Eusebio Valero.  A review on design of experiments and surrogate models in aircraft real-time and many-query aerodynamic analyses.  Progress in Aerospace Sciences, 96:23–61, January 2018.
  * [92] Dauda Ibrahim, Megan Jobson, Jie Li, and Gonzalo Guillén-Gosálbez.  Optimization-based design of crude oil distillation units using surrogate column models and a support vector machine.  Chemical Engineering Research and Design, 134:212–225, June 2018\.
  * [93] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  Nature, 521(7553):436–444, May 2015.
  * [94] Michael L. Thompson and Mark A. Kramer.  Modeling chemical processes using prior knowledge and neural networks.  AIChE Journal, 40(8):1328–1340, 1994.
  * [95] Chigozie Nwankpa, Winifred Ijomah, Anthony Gachagan, and Stephen Marshall.  Activation Functions: Comparison of trends in Practice and Research for Deep Learning, November 2018.
  * [96] D. G. Krige.  A statistical approach to some basic mine valuation problems on the Witwatersrand.  Journal of the Southern African Institute of Mining and Metallurgy, 52(6):119–139, December 1951.
  * [97] Carl Edward Rasmussen and Christopher K. I. Williams.  Gaussian processes for machine learning.  Adaptive computation and machine learning. MIT Press, Cambridge, Mass, 2006.
  * [98] Eddie Davis and Marianthi Ierapetritou.  A kriging method for the solution of nonlinear programs with black-box functions.  AIChE Journal, 53(8):2001–2012, aug 2007.
  * [99] Christopher M. Bishop.  Pattern recognition and machine learning.  Information science and statistics. Springer, New York, 2006.
  * [100] Alex Durkin.  Object-orientated derivative-free optimisation.  https://github.com/atdurkin/oodx, 2022.
  * [101] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [102] Geoffrey S. Simate, John Cluett, Sunny E. Iyuke, Evans T. Musapatika, Sehliselo Ndlovu, Lubinda F. Walubita, and Allex E. Alvarez.  The treatment of brewery wastewater for reuse: State of the art.  Desalination, 273(2):235–247, June 2011.
  * [103] Hashem Asgharnejad, Ehsan Khorshidi Nazloo, Maryam Madani Larijani, Nima Hajinajaf, and Hamidreza Rashidi.  Comprehensive review of water management and wastewater treatment in food processing industries in the framework of water-food-environment nexus.  Comprehensive Reviews in Food Science and Food Safety, 20(5):4779–4815, 2021.
  * [104] A. Gangagni Rao, T. Sasi Kanth Reddy, S. Surya Prakash, J. Vanajakshi, Johny Joseph, and P. N. Sarma.  pH regulation of alkaline wastewater with carbon dioxide: A case study of treatment of brewery wastewater in UASB reactor coupled with absorber.  Bioresource Technology, 98(11):2131–2136, August 2007.
  * [105] Mohd Najibul Hasan, Mohd Musheer Altaf, Nadeem A. Khan, Afzal Husain Khan, Abid Ali Khan, Sirajuddin Ahmed, P. Senthil Kumar, Mu. Naushad, Anushka Upamali Rajapaksha, Jibran Iqbal, Vineet Tirth, and Saiful Islam.  Recent technologies for nutrient removal and recovery from wastewaters: A review.  Chemosphere, 277:130328, August 2021.
  * [106] WRAP.  Compost and Digestate in Agriculture: Good Practice Guide, 2016\.
  * [107] Selami Demir.  Comparison of performances of biological nutrient removal systems for municipal wastewater treatment.  Sigma Journal of Engineering and Natural Sciences, 38(3):1235–1248, October 2021.  Number: 3.
  * [108] Jin Li, Chao He, Tian Tian, Zongkuan Liu, Zhaolin Gu, Guan Zhang, and Wendong Wang.  UASB-modified Bardenpho process for enhancing bio-treatment efficiency of leachate from a municipal solid waste incineration plant.  Waste Management, 102:97–105, February 2020.
  * [109] Maria Concetta Tomei, Valentina Stazi, Saba Daneshgar, and Andrea G. Capodaglio.  Holistic Approach to Phosphorus Recovery from Urban Wastewater: Enhanced Biological Removal Combined with Precipitation.  Sustainability, 12(2):575, January 2020.
  * [110] Elham Ashrafi, Arjomand Mehrabani Zeinabad, Seyed Mehdi Borghei, Elena Torresi, and Julian Muñoz Sierra.  Optimising nutrient removal of a hybrid five-stage Bardenpho and moving bed biofilm reactor process using response surface methodology.  Journal of Environmental Chemical Engineering, 7(1):102861, February 2019.
  * [111] Nikolaos V. Sahinidis.  BARON: A general purpose global optimization software package.  Journal of Global Optimization, 8(2):201–205, March 1996.
  * [112] Gonzalo Diaz, Achille Fokoue, Giacomo Nannicini, and Horst Samulowitz.  An effective algorithm for hyperparameter optimization of neural networks, May 2017.
  * [113] Diederik P. Kingma and Jimmy Ba.  Adam: A Method for Stochastic Optimization.  arXiv:1412.6980 [cs], January 2017.

## Supplementary information

### Neural network training

Training NNs consists of iterative algorithms to optimise the weights $W$, and
biases b for a given network structure. Training proceeds by making
predictions using the initialised value for the weights and bias parameters.
As such, these initial prediction are not very accurate. However, predictions
are made iteratively during the training cycle referred to as the forward
pass. One forward pass is completed by evaluating Equations 8a-8d for some
provided inputs. Each forward pass yields thereby yields predictions from the
output layer.

After each forward pass, the NN error is calculated by evaluating the loss
function. The loss function is used as the objective function in an
optimisation problem to minimise the error between NN predictions and training
observations. The loss function used depends on the data and problem that the
NN is being used to solve. Specifically, NNs for regression are trained using
loss functions such as the mean absolute error or mean squared error (Equation
42) to minimise the loss between $n$-vectors for network predictions
$\hat{\textbf{y}}$ and continuous regression observations y. Equation 43 shows
an NN cost function incorporating a loss function as well as an L2
regularisation term which sums the elements of the weight matrices, where
$\delta$ is a weight decay parameter and where the dependence on the weights
and biases $C(W,\textbf{b})$ is incurred via the predictions in the loss
function.

$L(\hat{\textbf{y}},\textbf{y})=\frac{1}{n}\sum_{i=1}^{n}\left(\hat{y}_{i}-y_{i}\right)^{2}$
(42)

$C(W,\textbf{b})=L(\hat{\textbf{y}},\textbf{y})+\frac{\delta}{2n}\sum_{i=1}^{n}w_{i}^{2}$
(43)

After calculating the error incurred in the forward pass, the backward pass
proceeds to calculate the changes in weights that provide the best reduction
in the cost function. This step is called backward propagation and typically
uses a gradient decent optimisation algorithm to determine the best direction
(changes in weights) to minimise the error. Backward propagation occurs in the
opposite direction to the forward pass, proceeding from the output later to
the input layer. The weights and biases are updated in the direction which
minimises the NN cost function which is determined by calculating the gradient
of the cost function with respect to the weights and biases. With respect to
the weights, the gradient of the cost function equals the gradient of the loss
function plus the gradient of the L2 weight regularisation term (Equation 44).
The gradient of the cost function is equal to the gradient of the loss
function with respect to the biases (Equation 45).

$\frac{\partial C}{\partial w^{(\lambda)}_{j,i}}={}\frac{\partial L}{\partial
w^{(\lambda)}_{j,i}}+\delta w^{(\lambda)}_{j,i}$ (44)

$\frac{\partial C}{\partial b^{(\lambda)}_{j}}={}\frac{\partial L}{\partial
b^{(\lambda)}_{j}}$ (45)

During back propagation, the gradients of the loss function with respect to
the weights can be determined by using the chain rule to propagation back
through the network (Equation 46).

$\displaystyle\frac{\partial L}{\partial w^{(\lambda)}_{j,i}}={}$
$\displaystyle\;\frac{\partial L}{\partial z^{(\lambda)}_{j}}\frac{\partial
z^{(\lambda)}_{j}}{\partial w^{(\lambda)}_{j,i}}$ (46) $\displaystyle={}$
$\displaystyle\;\frac{\partial L}{\partial a^{(\lambda)}_{j}}\frac{\partial
a^{(\lambda)}_{j}}{\partial z^{(\lambda)}_{j}}\frac{\partial
z^{(\lambda)}_{j}}{\partial w^{(\lambda)}_{j,i}}$

The derivative of the loss function with respect to a node output
$a^{(\lambda)}_{j}$ can be written as shown by Equation 47 where the
derivative of node input $z^{(\lambda+1)}_{k}$ with respect node output
$a^{(\lambda)}_{j}$ can be determined as $w^{(\lambda)}_{k,j}$ by
differentiating Equation 8a with respect to layer outputs.

$\displaystyle\frac{\partial L}{\partial a^{(\lambda)}_{j}}={}$
$\displaystyle\;\sum_{k=1}^{N_{\lambda+1}}\frac{\partial L}{\partial
z^{(\lambda+1)}_{k}}\frac{\partial z^{(\lambda+1)}_{k}}{\partial
a^{(\lambda)}_{j}}$ (47) $\displaystyle={}$
$\displaystyle\;\sum_{k=1}^{N_{\lambda+1}}\frac{\partial L}{\partial
z^{(\lambda+1)}_{k}}w^{(\lambda)}_{k,j}$

Equation 48 shows the derivative of node output $a^{(\lambda)}_{j}$ with
respect to node input $z^{(\lambda)}_{j}$ which can be determined as the
derivative of the activation function $\xi^{(\lambda)^{\prime}}$ by
differentiating Equation 8b with respect to layer inputs. Some derivatives of
common activation functions are shown in Table 1. Additionally, Equation 49
shows how the derivative of node input $z^{(\lambda)}_{j}$ with respect to
weight $w^{(\lambda)}_{j,i}$ can be written as node output
$a^{(\lambda-1)}_{i}$ by differentiating Equation 8a with respect to weights.

$\frac{\partial a^{(\lambda)}_{j}}{\partial
z^{(\lambda)}_{j}}=\xi^{(\lambda)^{\prime}}\left(z^{(\lambda)}_{j}\right)$
(48)

$\frac{\partial z^{(\lambda)}_{j}}{\partial
w^{(\lambda)}_{j,i}}=a^{(\lambda-1)}_{i}$ (49)

Finally, an error signal of node $j$ in layer $\lambda$,
$\epsilon_{j}^{(\lambda)}$ is defined by Equation 50 to simplify Equation 51
after substituting Equations 47-49 into Equation 46.

$\displaystyle\epsilon_{j}^{(\lambda)}={}$ $\displaystyle\;\frac{\partial
L}{\partial z^{(\lambda)}_{j}}$ (50) $\displaystyle={}$
$\displaystyle\;\frac{\partial L}{\partial a^{(\lambda)}_{j}}\frac{\partial
a^{(\lambda)}_{j}}{\partial z^{(\lambda)}_{j}}$ $\displaystyle={}$
$\displaystyle\;\sum_{k=1}^{N_{\lambda+1}}\left(\frac{\partial L}{\partial
z^{(\lambda+1)}_{k}}w^{(\lambda)}_{k,j}\right)\xi^{(\lambda)^{\prime}}\left(z^{(\lambda)}_{j}\right)$

$\displaystyle\frac{\partial L}{\partial w^{(\lambda)}_{j,i}}={}$
$\displaystyle\;\sum_{k=1}^{N_{\lambda+1}}\left(\frac{\partial L}{\partial
z^{(\lambda+1)}_{k}}w^{(\lambda)}_{k,j}\right)\xi^{(\lambda)^{\prime}}\left(z^{(\lambda)}_{j}\right)a^{(\lambda-1)}_{i}$
(51) $\displaystyle={}$
$\displaystyle\;\epsilon_{j}^{(\lambda)}a^{(\lambda-1)}_{i}$

The gradient of the cost/loss function with respect to the biases is similarly
approached by using the chain rule (Equation 52). The definition of the error
signal (Equation 50) as well as the derivative of Equation 8a with respect to
the biases (equal to 1), yields the simplified formulation of Equation 52

$\displaystyle\frac{\partial C}{\partial b^{(\lambda)}_{j}}={}$
$\displaystyle\;\frac{\partial L}{\partial z^{(\lambda)}_{j}}\frac{\partial
z^{(\lambda)}_{j}}{\partial b^{(\lambda)}_{j}}$ (52) $\displaystyle={}$
$\displaystyle\;\epsilon_{j}^{(\lambda)}$

Once the direction of steepest decent has been determined by an optimisation
algorithm, the NN parameters are updated by a given step size, $\alpha$, known
as the learning rate (Equations 53 and 54). The learning rate is another
hyperparameter which can be tuned by the user between the values 0 and 1: too
large and the optimiser might converge too quickly to a suboptimal local
solution; too small and the optimiser could get stuck and fail to find an
optimum within finite time. Optimisation algorithms available include
Stochastic Gradient Descent (SGD), L-BFGS, and Adam [113].

$\displaystyle w^{(\lambda)}_{j,i}={}$ $\displaystyle
w^{(\lambda)}_{j,i}-\alpha\frac{\partial C}{\partial w^{(\lambda)}_{j,i}}$
(53) $\displaystyle={}$
$\displaystyle\;w^{(\lambda)}_{j,i}-\alpha\left(\epsilon_{j}^{(\lambda)}a^{(\lambda-1)}_{i}+\delta
w^{(\lambda)}_{j,i}\right)$

$\displaystyle b^{(\lambda)}_{j}={}$
$\displaystyle\;b^{(\lambda)}_{j}-\alpha\frac{\partial C}{\partial
b^{(\lambda)}_{j}}$ (54) $\displaystyle={}$
$\displaystyle\;b^{(\lambda)}_{j}-\alpha\epsilon_{j}^{(\lambda)}$

Together, the forward and backward pass make one iteration which is typically
carried out for a subset of the training data called a mini-batch (where the
term batch is reserved for the whole data set). For a batch size equal to the
size of the training data set, this is referred to as batch gradient decent,
whereas for a batch size equal to 1 (i.e., 1 sample at a time) is what is used
specifically during SGD. For any batch size greater than 1 but less than the
size of the training set, the method is called mini-batch gradient decent.
Batch gradient decent is useful for fitting to relatively smooth functions
where the entire data can be used effectively to approximate a generalised
smooth output. However, batch gradient decent is computationally expensive for
very large data sets as gradients must be calculated for all samples at the
same time. On the other hand, SGD only calculates one gradient at a time and
takes an optimisation step accordingly, the result is quicker convergence for
very large data sets, but the global optimum can be difficult to find as the
objective function value fluctuates with every sample. Mini-batch gradient
decent achieves a trade-off between the two extreme methods, using small
batches enables better exploration of the solution space compared to SGD as
well as better computational times compared to batch gradient decent. Table 13
shows a pseudocode algorithm to train a NN by optimising the weights and
biases $W^{(\lambda)},\textbf{b}^{(\lambda)}$ via mini-batch gradient descent.

Table 13: Mini-batch gradient descent algorithm. input : $X$ (training
inputs), $Y$ (training observations), $n$ (number of training samples),
$\Lambda$ (number of hidden layers), $N^{(\lambda)}$ (number of nodes in each
layer), $\xi^{(\lambda)}$ (activation functions), epochs (number of epochs),
$n_{\text{batch}}$ (batch size), $\alpha$ (learning rate), $\delta$ (weight
decay), $L$ (loss function)
---
initialise : $W^{(\lambda)}\leftarrow W^{\text{init}}$ (random weights),
$\textbf{b}^{(\lambda)}\leftarrow\textbf{b}^{\text{init}}$ (random biases)
for epoch = 0 to epochs do |
Randomly shuffle training inputs $X$ | $\left.\begin{array}[]{l}\\\ \\\ \\\ \\\ \\\ \end{array}\right\\}\text{Mini-batch setup}$
Apply same permutation to training observations $Y$
for $i=0$ to $n$ step $n_{\text{batch}}$ do
$X^{\text{batch}}\leftarrow X[i:i+n_{\text{batch}}]$
$Y^{\text{batch}}\leftarrow Y[i:i+n_{\text{batch}}]$
Initialise empty array of predictions $\hat{Y}^{\text{batch}}$
for x in $X^{\text{batch}}$ do | $\left.\begin{array}[]{l}\\\ \\\ \\\ \\\ \\\ \\\ \\\ \end{array}\right\\}\text{Forward pass}$
$\textbf{z}^{(1)}\leftarrow W^{(0)}\textbf{x}+\textbf{b}^{(0)}$
for $\lambda=1$ to $\Lambda$ do
$\textbf{a}^{(\lambda)}\leftarrow\xi^{(\lambda)}\left(\textbf{z}^{(\lambda)}\right)$
$\textbf{z}^{(\lambda+1)}\leftarrow
W^{(\lambda)}\textbf{a}^{(\lambda)}+\textbf{b}^{(\lambda)}$
$\hat{\textbf{y}}\leftarrow\textbf{z}^{(\Lambda+1)}$
Update relevant row in $\hat{Y}^{\text{batch}}$ with $\hat{\textbf{y}}$
$C^{\text{batch}}\leftarrow L(\hat{Y},Y)+\frac{\delta}{2n}\sum_{i=1}^{n}w_{i}^{2}$ |  $\quad\;\triangleright$ Error evaluation
$\epsilon_{j}^{(\lambda)}\leftarrow\sum_{k=1}^{N_{\lambda+1}}\left(\frac{\partial L}{\partial z^{(\lambda+1)}_{k}}w^{(\lambda)}_{k,j}\right)\xi^{(\lambda)^{\prime}}\left(z^{(\lambda)}_{j}\right)$ |  $\quad\;\triangleright$ Backward propagation
for $\lambda=0$ to $\Lambda$ do | $\left.\begin{array}[]{l}\\\ \\\ \\\ \\\ \\\ \end{array}\right\\}\text{Optimisation step}$
for $k=0$ to $N_{\lambda+1}$ do
for $j=0$ to $N_{\lambda}$ do
$w_{k,j}^{(\lambda)}\leftarrow
w^{(\lambda)}_{j,i}-\alpha\left(\epsilon_{j}^{(\lambda)}a^{(\lambda-1)}_{i}+\delta
w^{(\lambda)}_{j,i}\right)$
$b^{(\lambda)}_{j}\leftarrow b^{(\lambda)}_{j}-\alpha\epsilon_{j}^{(\lambda)}$
return : $W^{(\lambda)},\textbf{b}^{(\lambda)}$ |

An epoch refers to a forward and backward pass of the entire data set (batch)
in multiple mini-batches, where an epoch contains a number of iterations equal
to the size of the batch divided by the size of the mini-batch. The number of
epochs should be tuned to ensure that the NN has enough exposure to the
training data to achieve a good fit, but not so much that it succumbs to
overfitting. To combat this, a validation data set can be used during training
to test the trained NN on unseen data (similar to the testing data set but
used during training to determine overfitting) and achieve regularisation.
During training, the training error decreases as the NN is updated and fitted
to the data. Initially, the validation error will also decrease as the NN goes
from initially underfit to being fitted. Eventually, as overfitting starts to
occur, the training error will continue to decrease but the validation error
will start to increase. This phenomenon can be used to stop training and avoid
overfitting, by specifying a large number of epochs along with a termination
condition based on the validation error increasing.

### OODX modelling objects

Table 14: DataHandler object. class DataHandler()
---
Modelling object to sample, process, and store data.
Attributes
space : array-like of shape (n_inputs, 2)
Lower and upper bounds on each input feature in the original space.
x : array-like of shape (n_samples, n_inputs)
Input samples in the original space.
y : array-like of shape (n_samples, n_outputs)
Output samples in the original space.
t : array-like of shape (n_samples, 1)
Binary convergence targets.
x_train : array-like of shape (n_train_samples, n_inputs)
Training input samples in the original space.
y_train : array-like of shape (n_train_samples, n_outputs)
Training output samples in the original space.
t_train : array-like of shape (n_train_samples, 1)
Training binary convergence targets.
x_test : array-like of shape (n_test_samples, n_inputs)
Testing input samples in the original space.
y_test : array-like of shape (n_test_samples, n_outputs)
Testing output samples in the original space.
t_test : array-like of shape (n_test_samples, 1)
Testing binary convergence targets.
space_ : array-like of shape (n_inputs, 2)
Lower and upper bounds on each input feature in the scaled modelling space.
x_ : array-like of shape (n_samples, n_inputs)
Input samples in the scaled space.
y_ : array-like of shape (n_samples, n_outputs)
Output samples in the scaled space.
x_train_ : array-like of shape (n_train_samples, n_inputs)
Training input samples in the scaled space.
y_train_ : array-like of shape (n_train_samples, n_outputs)
Training output samples in the scaled space.
x_test_ : array-like of shape (n_test_samples, n_inputs)
Testing input samples in the scaled space.
y_test_ : array-like of shape (n_test_samples, n_outputs)
Testing output samples in the scaled space.
x_mean : array-like of shape (n_inputs,)
Mean over n_samples of each input feature vector.
x_std : array-like of shape (n_inputs,)
Standard deviation over n_samples of each input feature vector.
y_mean : array-like of shape (n_outputs,)
Mean over n_samples of each output variable vector.
y_std : array-like of shape (n_outputs,)
Standard deviation over n_samples of each output variable vector.
x_train_mean : array-like of shape (n_inputs,)
Mean over n_train_samples of each input feature vector.
x_train_std : array-like of shape (n_inputs,)
Standard deviation over n_train_samples of each input feature vector.
y_train_mean : array-like of shape (n_outputs,)
Mean over n_train_samples of each output variable vector.
y_train_std : array-like of shape (n_outputs,)
Standard deviation over n_train_samples of each output variable vector.
Methods
init(n_samples, space, n_outputs=1, method=‘lhs’)
Generate initial input samples using static sampling strategies and initialise
empty arrays for outputs and binary convergence targets.
Parameters
n_samples : integer
Number of input samples to generate.
space : array-like of shape (n_inputs, 2)
Lower and upper bounds on each input feature in the original space.
n_outputs : integer, default=1
Number of output dimensions to initialise empty arrays.
method : string, default=‘lhs‘
Static sampling strategy to generate initial input samples.
Returns
self : object
DataHandler class instance.
split(test_size=0.3)
Train-test split on data.
Parameters
test_size : float, default=0.3
Fraction of samples to be assigned to testing sets.
Returns
self : object
DataHandler class instance.
scale()
Standardise input-output data. x is standardised and attributes the output as
x_. y is standardised to y_ using the mean and standard deviation of the
subset of feasible outputs to ensure that any infeasible outputs do not skew
the data. x_train is standardised to x_train_ using the mean and standard
deviation of the input training set. x_test is standardised to x_test_ using
the mean and standard deviation of the input training set to ensure that
testing data is mapped into the same space in which the model was trained.
y_train is standardised to y_train_ using the mean and standard deviation of
the subset of feasible training outputs to ensure that any infeasible outputs
do not skew the data. y_test is standardised to y_test_ using the mean and
standard deviation of the subset of feasible output training samples to ensure
that testing data is mapped into the same space in which the model was
trained. space is standardised to space_ using the mean and standard deviation
of the input training set is the data has been train-test split, otherwise
using the mean and standard deviation of the original input data.
Returns
self : object
DataHandler class instance.
scale_space(space)
Map new lower, upper bounds on the features into the scaled modelling space.
Parameters
space : array-like of shape (n_inputs, 2)
New lower and upper bounds on each input feature to be mapped into the scaled
modelling space.
Returns
new_space : array-like of shape (n_inputs, 2)
Lower, upper feature bounds mapped into the scaled modelling space.
scale_x(x)
Map new input samples into the scaled modelling space.
Parameters
x : array-like of shape (?, n_inputs)
New input samples to be mapped into the scaled modelling space.
Returns
scaled_x : array-like of shape (?, n_inputs)
New input samples mapped into the scaled modelling space.
inv_scale_x(x)
Map scaled inputs back into the original space.
Parameters
x : array-like of shape (?, n_inputs)
Input samples to be mapped back into the original space.
Returns
inv_scaled_x : array-like of shape (?, n_inputs)
Input samples mapped back into the original space.
scale_y(y)
Map new output samples into the scaled modelling space.
Parameters
y : array-like of shape (?, n_outputs)
New output samples to be mapped into the scaled modelling space.
Returns
scaled_y : array-like of shape (?, n_outputs)
New output samples mapped into the scaled modelling space.
inv_scale_y(y)
Map scaled outputs back into the original space.
Parameters
y : array-like of shape (?, n_outputs)
Output samples to be mapped back into the original space.
Returns
inv_scaled_y : array-like of shape (?, n_outputs)
Output samples mapped back into the original space.
Table 15: NN object. class NN(layers, activation=’tanh’, is_classifier=False)
---
Neural network modelling object for regression and classification.
Parameters
layers : list
List of the number of nodes in each layer.
activation : string, default=‘tanh’
Activation functions used throughout the network.
is_classifier : Boolean, default=False
Boolean variable equal to True if the NN is used for classification to change
the loss function to binary cross entropy loss on the logits.
Attributes
name : string
‘NN’ or ‘NNClf’ depending on whether is_classifier is False or True,
repectively.
layers : list
List of the number of nodes in each layer.
activation : string
Activation functions used throughout the network.
weights : list
List of weight matrices for each layer where the length of the list equals the
number of layers in the network ($\Lambda+1$). Each element of the list is
array-like of shape ($N_{\lambda}$, $N_{\lambda-1}$) where $N_{\lambda}$ is
the number of nodes in layer $\lambda$, for layers $\lambda=1,...,\Lambda+1$.
biases : list
List of bias vectors for each layer where the length of the list equals the
number of layers in the network ($\Lambda+1$). Each element of the list is
array-like of shape ($N_{\lambda}$,) for layers $\lambda=1,...,\Lambda+1$.
Methods
fit(x, y, batch_size=10, epochs=1000, learning_rate=1e-2, weight_decay=0.0,
loss_func=nn.MSELoss(), iprint=False)
Fit NN object to input-output data.
Parameters
x : array_like of shape (n_samples, n_inputs)
Input samples used to train the NN.
y : array_like of shape (n_samples, n_outputs)
Output samples used to train the NN.
batch_size : integer, default=10
Batch size used in batch gradient descent.
epochs : integer, default=1000
The number of iterations of forward and backward propagation.
learning_rate : float, default=0.01
The step size used during NN parameter optimisation.
weight_decay : float, default=0.0
Weight decay parameter used to regularise the NN during parameter
optimisation.
loss_func : object, default=torch.nn.MSELoss()
The loss function to be used during back propagation. Automatically switched
to binary cross entropy loss on the logits if is_classifier = True.
iprint : Boolean, default=False
If True, prints the time taken to train the NN.
Returns
self : object
Trained NN class instance.
predict(x)
Use a trained NN object to make predictions at new input data.
Parameters
x : array_like of shape (?, n_inputs)
New input samples used to make predictions.
Returns
predictions : array_like of shape (?, n_outputs)
Predictions from the NN made at new inputs.
formulation(x)
Use a trained NN object to make predictions at new input data based on an
explicit mathematical formulation.
Parameters
x : array_like of shape (?, n_inputs)
New input samples used to make predictions.
Returns
predictions : array_like of shape (?, n_outputs)
Predictions from the explicit mathematical formulation of the NN made at new
inputs.
Table 16: GPR object. class GPR(kernel=‘rbf’, noise=1e-10, porder=2)
---
Gaussian process regression model object.
Parameters
kernel : string, default=‘rbf’
Name of the kernel function to be used.
noise : float, default=1e-10’
Nugget term for modelling noisy observations.
porder : integer, default=2
Order of polynomial kernel function.
Attributes
name : string
Class attribute always equal to ‘GPR’.
kernel_name : string
Name of the kernel function implemented.
noise : float
Nugget parameter value.
x_train : array-like of shape (n_samples, n_inputs)
Training inputs as GPR model parameter.
length_scale : float
Optimised parameter for ‘rbf’ kernel.
constant_value : float
Optimised parameter for all kernel functions.
sigma_0 : float
Optimised parameter for linear/polynomial kernel functions.
inv_K : array-like of shape (n_samples, n_samples)
Inverse covariance matrix between training data.
porder : integer
Order of polynomial kernel function.
Methods
fit(x, y, iprint=False)
Fit GPR object to input-output data.
Parameters
x : array_like of shape (n_samples, n_inputs)
Input training samples.
y : array_like of shape (n_samples, 1)
Output training samples.
iprint : Boolean, default=False
If True, prints the time taken to train the model.
Returns
self : object
Trained GPR object.
predict(x, return_std=False, return_cov=False)
Use a trained GPR object to make predictions at new input data.
Parameters
x : array_like of shape (?, n_inputs)
New input samples used to make predictions.
return_std : Boolean, default=False
If True, returns the standard deviation of predictions.
return_cov : Boolean, default=False
If True, returns the covariance of predictions.
Returns
predictions : array_like of shape (?, 1)
Predictions from the GPR made at new inputs.
std : array_like of shape (?, 1), optional
Standard deviation in predictions from the GPR made at new inputs.
cov : array_like of shape (?, 1), optional
Covariance in predictions from the GPR made at new inputs.
formulation(x, return_std=False)
Use an explicit mathematical formulation of a trained GPR object to make
predictions at new input data.
Parameters
x : array_like of shape (?, n_inputs)
New input samples used to make predictions.
return_std : Boolean, default=False
If True, returns the standard deviation of predictions.
Returns
predictions : array_like of shape (?, 1)
Predictions from the GPR made at new inputs.
std : array_like of shape (?, 1), optional
Standard deviation in predictions from the GPR made at new inputs.
Table 17: GPC object. class GPC()
---
Gaussian process classification model object.
Attributes
name : string
Class attribute always equal to ‘GPC’.
x_train : array-like of shape (n_samples, n_inputs)
Training inputs as GPC model parameter.
t_train : array-like of shape (n_samples, n_inputs)
Training targets as GPC model parameter.
l : float
Optimised length scale parameter.
sigma_f : float
Optimised process variance parameter.
delta : array-like of shape (n_samples, 1)
Optimised GPC parameter.
inv_P : array-like of shape (n_samples, n_samples)
Optimised GPC parameter.
Methods
fit(x, t, iprint=False)
Fit GPC object to input-output data.
Parameters
x : array_like of shape (n_samples, n_inputs)
Input training samples.
t : array_like of shape (n_samples, 1)
Training targets.
iprint : Boolean, default=False
If True, prints the time taken to train the model.
Returns
self : object
Trained GPC object.
predict(x, return_std=False, return_class=False, threshold=0.5)
Use a trained GPC object to make probability predictions that t=1 at new input
data.
Parameters
x : array_like of shape (?, n_inputs)
New input samples used to make predictions.
return_std : Boolean, default=False
If True, returns the standard deviation of the latent variable at predictions.
return_class : Boolean, default=False
If True, returns class prediction as well as probability predictions.
threshold : float, default=0.5
Probability threshold for mapping probabilities onto class predictions.
Returns
predictions : array_like of shape (?, 1)
Probability predictions that t=1 at new inputs.
std : array_like of shape (?, 1), optional
Standard deviation of the latent function.
formulation(x)
Use an explicit mathematical formulation of a trained GPC object to make
predictions at new input data.
Parameters
x : array_like of shape (?, n_inputs)
New input samples used to make predictions.
Returns
predictions : array_like of shape (?, 1)
Probability predictions from the GPC made at new inputs.
Table 18: OODXBlock object. class OODXBlock(model)
---
Generate Pyomo Block formulations for machine learning model objects.
Parameters
model : object
Machine learning model object i.e., NN, GPR, GPC.
Attributes
model : object
Machine learning model object i.e., NN, GPR, GPC.
formulation : object
Pyomo Block formulation.
Methods
get_formulation(return_std=False)
Get Pyomo Block formulation for model attribute.
Parameters
return_std : Boolean, default=False
If True, returns a mathematical program to model the GPR standard deviation.
Returns
formulation : object
Pyomo Block formulation.
Table 19: AdaptiveSampler object. class AdaptiveSampler(space)
---
Object used to generate Pyomo-based adaptive sampling formulations.
Parameters
space : array-like of shape (n_inputs, 2)
Adaptive sampling search space.
Attributes
space : array-like of shape (n_inputs, 2)
Adaptive sampling search space.
delaunay : object
SciPy Delauanay triangulation object.
Methods
max_gp_std(model)
Generate Pyomo-based formulation to maximise the uncertainty of a GPR model
for model-based exploration.
Parameters
model : object
GPR model object.
Returns
m : object
Pyomo formulation.
max_triangle(x, include_vertices=False)
Generate Pyomo-based formulation to maximise the size of Delaunay triangulated
regions for data-driven exploration.
Parameters
x : array-like of shape (n_samples, n_inputs)
Input data on which to triangulate.
include_vertices : Boolean, default=False
If True, vertices of the search space will be included in the triangulation.
Returns
m : object
Pyomo formulation.
modified_expected_improvement(model, y, sense)
Generate Pyomo-based formulation to maximise the modified expected improvement
for a GPR model for model-based exploitation.
Parameters
model : object
GPR model object.
y : array-like of shape (n_samples, 1)
Output data from which to determine the current best sample.
sense : string
‘max’ or ‘min’ depending on whether to exploit around the current maximal or
minimal sample.
Returns
m : object
Pyomo formulation.
exploit_triangle(x, y, sense, include_vertices=False)
Generate Pyomo-based formulation to maximise the size of Delaunay triangulated
regions connected to the current best sample for data-driven exploitation.
Parameters
x : array-like of shape (n_samples, n_inputs)
Input data on which to triangulate.
y : array-like of shape (n_samples, 1)
Output data from which to determine the current best sample.
sense : string
‘max’ or ‘min’ depending on whether to exploit around the current maximal or
minimal sample.
include_vertices : Boolean, default=False
If True, vertices of the search space will be included in the triangulation.
Returns
m : object
Pyomo formulation.