Our general goal is to use NN modeling to solve nonlinear regression problems. In the considered SOFC application, this task corresponds to approximating scalar-valued multivariate functions that represent the electric power of the fuel cell in terms of all relevant influence factors. As already stated in the introduction, we restrict ourselves to the case of shallow NNs with a suitably chosen internal structure (cf. Rauh et al., 2021) and an optimized number of hidden layer neurons. This kind of network structuring is similarly exploited by Lutter et al. (2019), where NNs are employed to describe mechanical systems with the help of an approach based on the Lagrange formalism.

In contrast to our shallow NN approach, deep learning techniques are also employed in the current applied NN research. Here, the focus is on probabilistic regressions in combination with importance sampling, noise contrastive estimation, and maximum likelihood estimation (Gustafsson et al., 2020; Andersson et al., 2021; Gedon et al., 2021). In the area of probabilistic NN techniques, the estimation of probability density functions in a data-based context should be mentioned as a further approach that might find applications in control-oriented system modeling. Because of their limited applicability in the frame of control for fuel cell systems, other deep learning approaches, for example, for solving automatic (image) recognition tasks, are not discussed here. In our future research, we might consider extending the interval shallow NN modeling using the deep learning regression techniques that employ a feedforward or NARX structure similar to that considered in this paper. In such a way, the proposed approach could be extended not only to handle a control optimization in a bounded-error framework but also as an alternative to stochastic model predictive control techniques (González Querubín et al., 2020) with underlying deep learning system representations.

The general aim of this paper is thus to solve a function approximation problem (with possible memory effects represented by means of autoregressive NN structures) supposing that

• measurements of the system inputs are available at certain time instants,

Classical methods for point-valued NN parameterization make use of well-known optimization techniques such as the Levenberg-Marquardt algorithm, conjugate-gradient approaches, or Bayesian regularization (Dan Foresee and Hagan, 1997; Nocedal and Wright, 2006). These tools perform well if feedforward NNs with sufficiently large numbers of hidden layer neurons are considered. However, excessive overparameterizations should be avoided since this might lead to overfitting or high sensitivity against input disturbances which both cause poor generalizability to unknown inputs.

To prevent under- and overfitting of the data by the used NN, we employ the singular value decomposition approach theoretically derived in Kanjilal et al. (1993) and investigated in detail from an application (fuel cell) perspective in Rauh et al. (2021) and Rauh (2021) to find reasonable numbers of neurons in the hidden layer(s). In addition, this procedure is applicable to eliminate system inputs that do not contribute essentially to the input–output relations to be identified.

To solve the tasks listed above, we denote the vector of NN inputs by q k ∈ R m , where k ∈ {1, … , k _max} corresponds either to the sample number in a list of possible input values or to the data at the kth time instant if the system inputs possess a temporal relationship. Assuming a shallow NN with a single hidden layer, the vector y k ∈ R n of outputs can be computed by y k = N q k , W 1 , W 2 , b 1 , b 2 = W 2 ⋅ g W 1 ⋅ q k + b 1 + b 2 , (3) where g ⋅ : R L ↦ R L is the activation function in the hidden layer, L is the number of hidden layer neurons and the output layer is linear. Here, W 1 ∈ R L × m and W 2 ∈ R n × L are weighting matrices of suitable dimensions with the bias vectors b 1 ∈ R L and b 2 ∈ R n .

In the frame of training feedforward NNs with the aim of function approximation, the vector g x of activation functions is usually chosen either as the so-called ReLU function g x = g x 1 ⋮ g x L = max 0 , x 1 ⋮ max 0 , x L (4) or as a sigmoid-type function g x = σ x 1 ⋮ σ x L . (5)

A typical choice for the latter is the hyperbolic tangent σ x i = tanh x i , (6)

(implemented in Matlab by the function tansig). Throughout this paper, we use the representation σ x i = 2 1 + e − 2 x i − 1 (7) to avoid overestimation caused by a dependency of the numerator and denominator terms on the same variable in the alternative representation of the hyperbolic tangent function.

Remark 1

Both ReLU and sigmoid functions are capable of accurately representing input–output relationships in the frame of nonlinear function evaluation. ReLU representations are preferred in many current research works on regression problems in system identification because of their lower computational cost since efficient linear algebra routines and more simple gradient expressions can be applied during training phases. However, we make use of sigmoid representations in the hidden network layer in our paper because they automatically provide bounded outputs to input values with large magnitudes. This holds not only for point-valued NN realizations but also for the interval-valued counterparts introduced in this contribution. For those, constant interval diameters can be expected for inputs with large magnitude if sigmoid activation functions are used while the interval diameters would grow inevitably in the ReLU case because the vector components of (Eq. 4) are unbounded for x _i →+∞.

The NN model according to (Eq. 3) can be extended with the help of intervals if the inputs q _k are replaced by corresponding interval bounds to analyze the effect of input uncertainty on the network outputs. In this case, the NN parameters b ₁, b ₂, W ₁, and W ₂ remain fixed to the point values from the training phase. Related techniques for set-based reachability analysis were investigated in Tran et al. (2019a), Tran et al. (2019b), Xiang et al. (2018), Tjeng et al. (2019) and the references therein.

We demonstrate how this approach works using the following artificial example. Consider the task of approximating the static nonlinearity y = sin q 2 (8) by a feedforward NN. For that purpose, we generate N = 500 normally distributed input samples q _k with zero mean and standard deviation σ = 3 as the system inputs for which output measurements are available. These outputs y _m,k , where the subscript m stands for measured data, are corrupted by uniformly distributed, additive noise from the interval Δ ⋅ − 1 ; 1 with Δ = 0.15. The presence of such noise processes is quite common in practice. It models the effects of random disturbances that can be considered as bounded due to quantization effects and the use of intelligent sensors detecting measurement outliers ¹ .

Figure 1 shows a comparison of the exact system behavior (Eq. 8), the training samples, and the corresponding point-valued NN representation, where the number of hidden layer neurons is set to L = 10 if either sigmoid-type or ReLU activation functions are used. In both cases, the NN was trained using a Bayesian regularization technique from Matlab. During this process, the overall data set was subdivided randomly into training (70%), test (15%), and validation (15%) data.

A naive interval extension of the network (i.e., that directly replacing all floating point operations with the corresponding interval counterparts) for an equidistant interval mesh q ̲ i ; q ̄ i of NN inputs with either q ̄ i − q ̲ i = 0.1 or q ̄ i − q ̲ i = 0.01 , where min ( q ̲ i ) = − 10 and max ( q ̄ i ) = 10 is shown in Figure 2. This interval extension aims at quantifying the ranges of possible NN outputs without the necessity for sampling the inputs in a point-valued form.

As mentioned in Remark 1, the use of a sigmoid hidden layer ensures bounded outputs for inputs with large absolute values. For both kinds of activation functions, the comparison of different mesh sizes for the input intervals reveals the problem of multiple interval dependencies that leads to overestimation in the computed output ranges. For that reason, it is necessary to apply techniques for the reduction of overestimation during the interval evaluation of NNs with interval-valued input data q _k . This might include a (brute force) subdivision of input intervals ² as shown in Figure 2 or algorithmic means such as centered form or affine representations and slope arithmetic (Cornelius and Lohner, 1984; Chapoutot, 2010). Due to the restriction to globally defined (and, for the sigmoid case, differentiable) functions in the hidden layer, these methods can be implemented using existing interval libraries such as IntLab (Rump, 1999) or JuliaIntervals (Ferranti, 2021). An advantage of the latter is that it can be employed to execute the NN’s source code directly on a GPU (graphics processing unit) to exploit data parallelism for large numbers of input samples q _k .

In contrast to performing a set-based reachability analysis of NN models that are trained with the help of local optimization procedures and thus possess point-valued parameters b ₁, b ₂, W ₁, and W ₂ as summarized in the previous subsection, the main goal of our paper is to determine suitable interval versions of these parameters to quantify the NNs’ modeling quality.

Gradient-based training techniques can be generalized quite easily to determine interval-valued NN parameterizations (Oala et al., 2020). However, a drawback of the direct NN training using interval parameters might be the convergence to poor local minima of the considered cost function if inflating layer weights leads to similar output intervals as replacing selected bias terms by interval parameters. During the search for interval-valued NN parameterizations, the cost function should not only reproduce the training data with high accuracy but also reduce the interval diameters of the network’s outputs, see Sec. 3 for further details. Additionally, we have to make sure that all (point-valued) system outputs (i.e., the actual point-valued measured data available for the training) are fully included in the networks’ output intervals. This property can be relaxed, for example, by allowing a certain percentage of outliers that do not necessarily need to be included in the network’s output range. This modification is similar to that from the relaxed intersection approach described in Jaulin and Bazeille (2009).

Poor convergence of the interval-based NN training relying on local optimization techniques by a direct search for interval parameters as stated above is caused by the fact that the conversion of different network parameters to interval quantities might lead to ambiguous output behavior, at least for some of the domains of the inputs q _k . This fact is visualized in Figure 3 for the NN from the previous subsection. To limit the overestimation caused by interval-valued inputs q k ∈ q k , the evaluation is performed on the fine grid from Figure 2B. To illustrate the influence of converting either b ₁, b ₂, W ₁, or W ₂ into interval quantities, each of them is multiplied independently with the interval bounds 0.95 ; 1.05 (parameter uncertainty of ±5%). Since b ₁, b ₂, W ₁, and W ₂ each appear only once in Eq. 3, the NN does not exhibit classical dependency-based overestimation. However, the wrapping effect described in Jaulin et al. (2001) might still be present where the matrix multiplication is involved. Output ambiguities, concerning a similar shape of the enclosures y k , can be seen when comparing Figure 3A with Figure 3C (and Figure 3B with Figure 3D). This visual comparison indicates that converting all entries in W _i to intervals might have a similar influence as the interval conversion of b _i for some of the respective entries.

To circumvent such kind of poor convergence and to reduce the effect of the illustrated ambiguities, we propose a two-stage identification procedure for feedforward NNs with interval parameters. Rather than just replacing all entries by intervals provided by local optimization in a single stage, we show in Section 3 how to specifically determine those entries in b _i and W _i in a two-stage procedure for which the interval conversion leads to the smallest diameters of the NN output intervals while nonetheless including (all) measured data.

In contrast to static function approximations, dynamic NN models contain a kind of memory for previous input and state information. This can be realized, for example, by training recurrent networks with an Elman or Jordan structure (Elman, 1990; Mandic and Chambers, 2001), training NARX models (Nelles, 2020), or combining static function approximation networks with dynamic elements that can be represented by integer- or fractional-order transfer functions. The latter option leads to the dynamic NN models and Hammerstein-type nonlinearity representations discussed in Rauh (2021).

In all cases, interval parameters within the corresponding networks or also interval-valued system inputs and initial conditions lead to discrete- or continuous-time dynamic systems with interval parameters. For such systems, the wrapping effect of IA is an important issue to handle. For corresponding approaches, the reader is referred to the vast literature on verified simulation procedures for dynamic systems, cf. Nedialkov (2011); Lohner (2001). As far as dynamic NN-based model evaluations are concerned, Xiang et al. (2018) and the references therein provide a good overview of the current state-of-the-art. However, the computational effort increases by a large degree in comparison to the evaluation of purely static NN models.

For that reason, we do not make use of a direct IA evaluation of any of the mentioned dynamic NN models. Instead, a NARX model, representing a point-valued approximation of the dynamic system behavior, is combined with a static interval-valued additive error correction model, both relying on the same finite window of input and output information.

For the purpose of static function approximations, an interval NN model structure shown in Figure 4 is considered. A generalization of the representation given in Eq. 3, namely, its interval extension aiming at a guaranteed enclosure of all measured samples y _m,k , is denoted by y m , k ∈ N q k , W 1 , W 2 , b 1 , b 2 = W 2 ⋅ g W 1 ⋅ q k + b 1 + b 2 , (9) where the interval parameters, highlighted by the thick lines in Figure 4, are defined for i ∈ {1, 2} according to b i = b ̌ i + [ Δ b ̲ i ; Δ b — i ] and W i = W ̌ i + [ Δ W ̲ i ; Δ W W — i ] . (10)

In analogy to the definition of scalar interval variables, the inequalities Δ b ̲ i ≤ Δ b ̄ i and Δ W ̲ i ≤ Δ W W — i hold in an element-wise manner. To make it possible to combine a point-valued NN parameterization with an IA-based uncertainty model, we assume that both b ̌ i and W ̌ i denote the interval midpoints to be identified in the first stage of the following procedure, while the additive bounds in the relations stated in (Eq. 10) are determined in the second stage.

In Figure 4 and all subsequent block diagrams of NN models in this paper, the parameters indicated by solid arrow heads are determined during the training and optimization phases of the proposed procedure while the non-filled arrow heads denote connections with fixed weights that are set to the constant value one.

We suggest the following interval procedure for static NN modeling.

Stage A Train the NN N q k , W ̌ 1 , W ̌ 2 , b ̌ 1 , b ̌ 2 according to (Eq. 3) by means of a standard algorithm (e.g., trainbr in Matlab) to minimize the quadratic cost function J N = ∑ k ∈ T ∑ i = 1 n y i , k − y m , i , k 2 , (11) where T denotes the set of training samples’ indices. The resulting parameters b ̌ 1 , b ̌ 2 , W ̌ 1 , W ̌ 2 define the interval midpoints that are kept constant during the following stage.

Stage B Determine the NN’s interval extension.

B1 Compute interval correction bounds for the parameter b ̌ 2 according to Δ b ̲ 2 ; Δ b ̄ 2 = y m , k − ⨆ k = 1 k max N q k , W ̌ 1 , W ̌ 2 , b ̌ 1 , b ̌ 2 (12) so that all measured ouput samples y _m,k are included in the interval-valued NN outputs. In (Eq. 12), the symbol ⨆ denotes the tightest axis-aligned interval enclosure around all arguments of this operator.

B2 Determine the intermediate minimum cost function value J * = ∑ k = 1 k max 1 T ⋅ w ( N ( q k , W ̌ 1 , W ̌ 2 , b ̌ 1 , b ̌ 2 + [ Δ b ̲ 2 ; Δ b — 2 ] ) ) (13) with w x being the element-wise extension of the interval diameter definition from (Eq. 1) and the vector 1 = 1 … 1 T ∈ R n . (14) B3 Minimize the cost function J = ∑ k = 1 k max 1 T ⋅ w N q k , W 1 , W 2 , b 1 , b 2 + J * ⋅ P k (15) by searching for optimal parameterizations of the additive interval bounds Δ b ̲ i ; Δ b ̄ i and Δ W ̲ i ; Δ W ̄ i in (Eq. 10). Here, P _k is a sufficiently large penalty term P k = 100 ⋅ δ k + 2 2 + δ k − 2 2 + Δ y k + 2 2 + Δ y k − 2 2 (16) with the squared Euclidean norms x 2 2 preventing minima that are inadmissible due to a violation of the network’s output intervals by individual measurements. These points are detected with the help of the Boolean indicator vectors δ k + = y m , k − y — k > 0 and δ k − = y m , k − y ̲ k < 0 (17) with their corresponding element-wise defined excess widths Δ y k + = max ( 0 , y m , k − y — k ) and Δ y k − = min 0 , y m , k − y ̲ k . (18)

Remark 2

Using this minimization procedure, implemented by means of a particle swarm optimization technique (Kennedy and Eberhart, 1995; Yang, 2014), it is guaranteed that the optimized bounds in Step B3 become tighter than the heuristic outer bounds in Step B1 in those domains for the inputs q _k, where dense identification samples are available. This is a direct consequence of the fact that the relation J ≤ J* between the cost functions (Eq. 13) and (Eq. 15) is always satisfied if P _k = 0. In this case, the NN is parameterized so that none of the measured samples lies outside the output interval bounds.

Remark 3

For scenarios in which a certain percentage of outliers are tolerated, the penalty term J* ⋅ P _k in (Eq. 15) can be modified easily to ignore those samples in the summation that have the largest indicator values P _k . For an example of this relaxation of the inclusion property, Sec. 4.

In Figure 5A and Figure 5B, a comparison is shown between the heuristic interval parameterization of a feedforward NN according to Step B1 and the systematic version according to Step B3 (for sigmoid activation functions). The comparison for the ReLU case is illustrated in Figures 5C,D.

These two scenarios make use of the same data as in the previous section. Due to the symmetric (uniformly distributed) nature of the disturbance, the interval parameterization can be simplified in this example by setting Δ b ̲ i = − Δ b ̄ i and Δ W ̲ i = − Δ W ̄ i (19) to reduce the degrees of freedom during the particle swarm optimization in Stage B. For both kinds of NN, it can be seen clearly that the interval diameters are much smaller in the domains in which dense measurements are available. Outside those domains, the sigmoid version of the NN has the advantage of providing tighter bounds due to the inherited saturation of the hyperbolic tangent function (Eq. 7). For that reason, we use only the sigmoid implementation in the context of SOFC modeling (cf. Sec. 4).

In Table 1, there is a summary of the computing times ³ for evaluating the trained neural networks. Obviously, the slowest version is the naive evaluation of the point-valued networks in the form of the Matlab network structure using the command sim. This can be accelerated significantly by employing automatic code generation (genFunction with the option MatrixOnly), which, however, pre-evaluates all NN-internal relations with the point-valued parameters. To obtain an implementation that is directly suitable for replacing the NN parameters with intervals, an extension of the Matlab functionalities by a vectorized implementation is necessary to which all parameters are passed as function arguments. Compared with this version, the interval-valued realization increases computing times by a factor 5 for ReLU activation functions and by a factor 10 in the nonlinear, sigmoid case.

Remark 4

Note that the true function (usually not available for the identification procedure) is not necessarily fully included in the determined interval outputs if the number of training points is insufficient in certain domains. If this phenomenon is detected at runtime after new measured data arrive, the model can be adapted by a re-application of Step B1. If computational resources allow this, Step B3, hot-started with the previously identified interval bounds, can also be re-applied.

In this paper, we extend the two different system structures for NARX models depicted in Figure 6 using IA. The first is the dynamic NN model N d,A in which both input variables q _k and output information from a (previous) time window influence the dynamics in a nonlinear way (Figure 6A). In contrast, the network N d,B in Figure 6B contains a simplified version in which only the input variables enter the nonlinear hidden layer with the sigmoid activation functions (Eqs. 5–7), while the previous outputs are fed back in a linear manner. For a compact mathematical representation of the networks y k = N d,A q k − M : k , y k − M : k − 1 , W d , q , 1 , W d , y , 1 , W d , 2 , b d , 1 , b d , 2 = W d , 2 ⋅ g W d , q , 1 ⋅ q k − M : k + W d , y , 1 ⋅ y k − M : k − 1 + b d , 1 + b d , 2 (20) and y k = N d,B q k − M : k , y k − M : k − 1 , W d , q , 1 , W d , y , 1 , W d , 2 , b d , 1 , b d , 2 = W d , 2 ⋅ g W d , q , 1 ⋅ q k − M : k + b d , 1 + W d , y , 1 ⋅ y k − M : k − 1 + b d , 2 (21) with M ≥ 1 as the number of previous sampling points, we introduce the stacked vectors q k − M : k = q k − M T … q k T T ∈ R m ⋅ M + 1 and y k − M : k − 1 = y k − M T … y k − 1 T T ∈ R n ⋅ M (22) for input and output variables, respectively.

As explained in Sec. 2.4, we suggest to extend both types of networks (ι ∈ {A, B} in Figure 6) using an additive correction that is implemented by means of a static NN with interval parameters such that the enclosure property y m , k ∈ N d , ι q k − M : k , y k − M : k − 1 , W d , q , 1 , W d , y , 1 , W d , 2 , b d , 1 , b d , 2 + N q k ′ , W 1 , W 2 , b 1 , b 2 , ι ∈ A , B (23) with the combined input vector q k ′ = q k − M : k T y k − M : k − 1 T T ∈ R m ⋅ M + 1 + n ⋅ M (24) holds for all measured samples y _m,k . Due to this specific structure, the network N with interval parameters can be optimized with the help of the two-stage procedure described in Sec. 3.1 as soon as the structure and the parameters of the point-valued NARX models N d , ι have been fixed.

A detailed comparison of the alternatives from Figure 6 is given in the following section using the example of modeling the electrochemical behavior of an SOFC system.

As the fundamental, static electric power model, the NN representation in Figure 7 is employed. This network is trained with the measured data described above by using the Matlab NN toolbox. We rely on the Bayesian regularization back-propagation algorithm (parallelized on four CPU cores) with a maximum number of 5,000 epochs to solve the training task. During the network training, a worsening of the validation performance has been allowed in 50 subsequent iterations with a random subdivision into training (70%), test (15%), and validation (15%) data. For further details, see Rauh (2021).

The visualization of the training results in Figure 9 reveals that the largest approximation errors occur during those operating phases of the SOFC stack in which either the hydrogen mass flow or the electric stack current show significant temporal variations. Therefore, point-valued NARX models (not yet considered in previous publications by the authors) are discussed in the following subsection to reduce approximation errors in dynamic operating phases.

As indicated in Figure 6, two different types of NARX models can be considered. These are either models with fully nonlinear characteristics of both the system inputs and the fed-back stack voltage ( N d,A ) or a model simplification in which only nonlinearities with respect to the input variables are accounted for ( N d,B ) . To ensure comparability between the previous static system model and the NARX representation, L = 7 is again used as the fixed number of hidden layer neurons. In addition, M = 10 (for the data set with the reduced sampling frequency of 1 Hz) serves as the memory length of the autoregressive model. This heuristically chosen value is motivated by a delay of a few seconds between the fuel cell stack’s inputs and outputs that became visible during the experimental data acquisition.

The training of both point-valued models N d,A (Figure 10) and N d,B (Figure 11) is based on a two stage procedure (using the same optimization algorithm as for the static case in the previous subsection). First, an open-loop training is carried out. Second, the NARX model obtained from the first stage is converted into a closed-loop structure by feeding back the simulated outputs. This configuration is then re-optimized to further reduce approximation errors and to avoid instability of the closed-loop NN models. Here, directions for our future research will be to optimize the memory length M and to develop techniques for automatically selecting different memory lengths for the network inputs q _k−M:k and for the output feedback y _k−M:k−1. When fully neglecting the output feedback, this would lead to nonlinear finite impulse response models in which the only delayed variables are the input data q _k−M:k .

Compared with the static NN model, both NARX options reduce the deviations between measured and simulated data significantly, especially during the dynamic operating phases. This is confirmed in Table 2 that provides information on the roots of the corresponding mean square approximation errors (RMS). It can be clearly seen that approximation errors for both the voltage and the power are approximately 30% lower if the dynamic model N d,B with the linear output feedback is used instead of the static representation. These errors can be reduced even further if the fully nonlinear representation N d,A is used instead (cf. the visualization of the approximation errors in Figures 10B, 11B).

In this section, we employ the procedure described in Sec. 3 to parameterize the interval extension of the static NN model. The left column of Figure 12 shows the results of the direct application of this procedure to the complete available data set for the SOFC system. It can be seen that the operating domains with the largest uncertainty of the model (characterized by the largest interval diameters of the forecast electric power in Figure 12C) are those regions in which significant power variations due to sharp changes of the fuel cell current or the hydrogen mass flow occur (cf. Figure 8). Here, symmetric bounds have been determined for all NN parameters as suggested in (Eq. 19).

To reduce the resulting interval widths for the static approximation of the electric fuel cell power, the penalty term (Eq. 16) is modified in such a way that the largest positive deviations in (Eq. 18) of y m , k − y ̄ k and y m , k − y ̲ k —corresponding to outliers—are ignored in both cases for at most 2% of the measured voltage samples (cf. the right column of Figure 12). This modification results in an interval extension of the NN model that encloses almost all measured voltage and power data, respectively. Note that outliers occur only at those points where step-like current changes of large amplitudes take place. If the resulting interval-valued model is employed in future work to improve the robustness of the maximum power point tracking procedures or online optimization approaches for enhancing the fuel efficiency according to Rauh et al. (2020) and Rauh (2021), such large current changes can be avoided easily by imposing suitable rate constraints on the control signals.

In the final part of this section, interval extensions of both possible NARX models are presented in Figure 13. As in the previous subsection, we again compare the cases with and without admissible outliers. It can be seen that the fully nonlinear NARX model N d,A does not only provide smaller approximation errors for the point-valued system identification but also yields significantly smaller diameters in its interval extension.

In contrast to the interval extension of the static NN model, it can be seen that the interval diameters in Figures 13E,F are almost proportional to the electric current. This indicates that the identified interval network parameters in (Eq. 23) lead to almost constant interval widths if the fuel cell stack voltage is concerned as the output of the dynamic system model. On the one hand, this confirms quite strongly that the point-valued system model is parameterized in terms of a fairly robust minimum of the employed quadratic error functional in Stage A of the optimization algorithm. On the other hand, this indicates that further improvements of the modeling accuracy could be achieved only by replacing the feedback of the fuel cell voltage in the NARX models by the electric power. However, this option is not further investigated in this paper, because it would require removing the physically motivated multiplicative system output shown in Figure 7. As observed in Rauh et al. (2021), this would lead to the necessity of a larger number of hidden layer neurons already during the phase of optimizing the point-valued system model. Therefore, a promising compromise for future research would be to modify the interval model (Eq. 23) in such a way that the point-valued NARX models N d , ι provide the stack voltage as a system output that is multiplied with the measured current I _k , while N ′ represents interval bounds for the power approximation error according to P el, m , k ∈ N d , ι q k − M : k , y k − M : k − 1 , W d , q , 1 , W d , y , 1 , W d , 2 , b d , 1 , b d , 2 ⋅ I k + N ′ q k ′ , W 1 , W 2 , b 1 , b 2 , ι ∈ A , B . (25)

It is especially interesting to investigate this modified formulation in an ongoing research, in which the Kalman filter-based online power identification according to Rauh (2021) and the interval-valued NN models from this paper are compared with respect to their capabilities to design robust procedures for the optimization of fuel efficiency while simultaneously tracking desired temporal responses for the stack power.