Methods

: The paper introduces a subclass of nonlinear differential-algebraic models of interest for applications. By restricting the nonlinearities to multilinear polynomials, it is possible to use modern tensor methods. This opens the door to new approximation and complexity reduction methods for large scale systems with relevant nonlinear behavior. The modeling procedures including composition, decomposition, normalization, and multi-linearization steps are shown by an example of a local energy system with a nonlinear electrolyzer, a linear buck converter and a PI controller with saturation.


Motivation
Automation systems in general are hybrid, i. e., some signals are discrete-valued and some are continuousvalued [20]. The discrete-valued part can be modeled by binary functions, which are generically in the multilinear subclass of polynomials. But multilinear functions can also be defined over real numbers, thus including linear functions, see Figure 1.
Linear modeling will not be sufficient for applications where larger deviations from an operating point have to be considered and which show essential nonlinear effects within their region of operation. Multilinear functions generically enable the usage of tensor algorithms and multilinear algebra is the mathematical domain of interest. In the last decade this domain has delivered results in the field of tensor decompositions which have been enabling break-throughs in many disciplines like quantum information theory or thermodynamic models of atomic alloy structures [5,25]. The research community of control and automation is only slightly aware of this: an open invited session was organized at the IFAC World Congress 2017, but e. g., the number of papers at IFAC WC 2020 dealing with tensors methods still was low. This paper follows the idea of using tensor decomposition techniques for parameter spaces of models. It further extends the class of multilinear time-invariant (MTI) models from explicit state space models to implicit descriptor models. This opens the door to reduced efficient modeling of huge nonlinear complex systems -where many challenging problems of today and in future are coming from.

Known results
Although tensor decompositions were defined a century ago [11] there are still open problems [3], some of them related to the field of algebraic geometry [6]. But applied mathematics and informatics research have found efficient strategies for approximative solutions [13]. The attention in engineering sciences has been growing with the availability of computational toolboxes [30,2,1,24]. The connection of tensor decomposition and multilinear systems is established [27].
Remark: Multiple or piecewise linear systems should be called multi-linear to distinguish them from multilinear systems. The first class needs switches and linear algebra, the latter involves multilinear algebra. Tensor representations can be found for both, but this paper is focussed on the second.
In the past years, explicit multilinear time-invariant (eMTI) models have been developed and used for control of HVAC systems in buildings [14,18,19,15]. The class of eMTI models is not closed w. r. t. basic compositions, but only the larger class of polynomial models [26]. Descriptor models which combine algebraic and differential equations are intrinsically composable. They describe the system dynamics implicitly and as well have been proven to be useful in many applications [22]. To the knowledge of the authors, extensions of eMTI models to implicit multilinear time-invariant (iMTI) models haven't been investigated so far. This will be the main and methodological contribution of this paper which will be illustrated by an example of an energy system.
Current and future problems of energy systems are the integration of new components in power networks -in the real world as well as in the models. Basic relevant literature for the energy system components discussed later can be found in [21]. The power system is undergoing a deep transformation towards power electronic dominated generation and consumption, which brings challenges regarding the safe operation of the grid. Therefore new techniques are needed to access the safe operation. The MTI class shows several aspects which makes it promising for power systems description, e. g., switching and discrete dynamics, multilinearity and scale.

Open questions
The non-closedness of eMTI models w. r. t. compositions has several drawbacks for complex system modeling which this paper overcomes by answering the questions: 1. How are basic connections (parallel, series, feedback) best represented? 2. Which model class holds the pros of multilinearity but offers composability? 3. Can complex energy system components be modeled within this class? 4. What complexity and accuracy do tensor energy network models have?
The paper is structured as follows: the next Section 2 recaps the current state of research of explicit multilinear modeling and introduces implicit models in matrix representation. Section 3 shows by an counterexample why explicit multilinear models are not closed w. r. t. composition of subsystems and gives the important result that implicit multilinear models are closed w. r. t. series, parallel and feedback compositions. In Section 4, tensor representations are given for explicit as well as implicit multilinear models together with one relevant decomposition method. How the canonical polyadic (CP) decomposed tensor representation can be further reduced by normalization is shown in Section 5. All methods are applied in Section 6 to a small scale energy system example for demonstration of the capabilities of the approach. The paper closes with conclusion and outlook.

Multilinear models
In this section the description of multilinear models is given. This class of models extends the linear class.

Multilinear functions
A function f (x) with x = (x 1 , . . . , x n ) is called multilinear, if the function is linear in each individual variable x i , meaning all other variables are held constant. This means, that multiplicative combinations of the scalars in x are feasible for the function to be multilinear. All combinations of the scalars can be generated by the so called monomial vector m(x), which is given by where ⊗ denotes the Kronecker product.
Example 2.1. For x ∈ ℝ 3 the monomial is given by Thus, the general form of a multilinear function can be given by the inner product i. e., the multiplication of a row vector f T of coefficients times the column vector m(x) of monomials. A very interesting property of multilinear functions is, that all Boolean functions belong to that class. In the following, the set denoted by = {0, 1} ≡ {FALSE, TRUE} refers to binary variables. Any Boolean function b of n variables can be represented by 2 n rows of a truth table, corresponding to n-element binary vectors given in lexicographical order. Assuming this order, the Boolean function can be represented by the truth vector b ∈ 2 n of the last column of the truth table. To show this, a vector of literals is introduced as the Kronecker product of the literals of the variables of the Boolean function, given as of order n is a multilinear square-free polynomial given by the scalar product of a Boolean truth vector b ∈ 2 n and a vector l(x) ∈ ℝ 2 n of literals. Zhegalkin Polynomials are sometimes also referred to as Reed-Muller Forms or Algebraic Normal Forms (ANF).
The following example of order n = 2, taken from [28], illustrates the idea.

Example 2.2.
For the Boolean function b = ¬(x 1 ∧ x 2 ) the truth table is given as As Zhegalkin polynomials (5) are multilinear, they can be equivalently represented by the parameter vector f of (3).
If Boolean values are inserted for x 1 and x 2 this results in Boolean values which correspond to the truth table.

Explicit multilinear models
With the monomial vector, an eMTI state space model can be written in continuous timė where u ∈ ℝ m is the input vector, y ∈ ℝ p the output vector, F ∈ ℝ n×2 (n+m) the transition matrix and G ∈ ℝ p×2 (n+m) the output matrix. The monomial vector of more than one vector treats the input vectors as concatenated, for two vectors i. e., m(x, u) = m((x 1 , x 2 , u 1 , u 2 )).

Remark.
In contrast to multilinear models, linear models are restricted to the subclass of linear functions going through the origin (f (0) = 0), whereas the first element of the monomial vector contains a constant.
Remark. Bilinear models are also a subclass of multilinear models where all monomials are restricted to order 2 and combinations of inputs and states.
Boolean functions can be represented as given in (5). If the Boolean variable is depended on a continuous variable x, the result of (5) will be continuous. Thus, continuous multilinear function and Boolean functions can be modeled individually by (6) and (7). On the other side, if continuous and Boolean functions are connected with each other a quantizer is needed, see [27] for more details. A simple solution is to use a Heaviside function as quantizer. Using this function, a saturation can be modeled as given in the following example.

Example 2.3.
Assuming the input to the saturation is a state, and the output is saturated between two values, e. g., the output of a controller needs to be between 0 and 1. The saturation function is depicted in Figure 2.
To model this as an eMTI model two discrete-valued variables are introduced The output is then calculated by with the slope s =

Implicit multilinear models
The class of explicit multilinear models has the drawback, that by connecting multiple eMTI models the overall model is not necessarily an eMTI model which will be shown in section 3. To overcome this obstacle we introduce iMTI models similar to descriptor systems in the next subsection.
For iMTI models a kernel representation is chosen and the monomial vector is extended. An iMTI model in kernel representation is given as where the states areẋ ∈ ℝ n , x ∈ ℝ n , the inputs are u ∈ ℝ m , the outputs are y ∈ ℝ p and the model matrix is H ∈ ℝ N×2 2n+m+p . The dimension of the model matrix is not only determined by the number of states, inputs and outputs but by the application specific number N of equations. A designated output equation is in implicit form not required, but output equations can be included in the model matrix H.

Remark.
Note that all eMTI models can be written as iMTI models by subtracting the right hand sides of (6) and (7) on both sides, leading to kernel representationṡ given as a semiexplicit differential algebraic equation, [4]. These equations can then be written in the form (12), because their left hand sides are obviously multilinear functions in m(ẋ , x, u, y).
Since the equations for the states are no longer in explicit form, the mapping from continuous-valued variables to discrete-valued variables can no longer be performed by quantization (8), but discrete valued states have to be defined in pairs, because the single kernel equation b(1−b) = b − b 2 = 0 of a Boolean constraint is polynomial, but not multilinear.

Example 2.4.
Assume the states b 1 and b 2 are Boolean and can therefore only take the values 0 or 1, meaning FALSE or TRUE, respectively. Introducing the equations forces at least one of the variables to be 0 by (15), but prevents both to be 0, because of (16), and thus forcing one to be 0 and the other one to be 1.
As the example shows, an iMTI model is capable to define Boolean variables. But to, e. g., model saturation, inequalities are no longer avoidable -which leads to the generic multilinear model Note that (17) includes the kernel representation (12) by setting and thus, the iMTI model (17) has the broadest spectrum of applicability. Next, saturation is modeled generically as iMTI.

Example 2.5.
To model saturation, we use the discrete variables b 1 and b 2 of Example 2.4 and append a second pair of discrete variables by Then the variables need to be connected to the ranges in Figure 2 given as In case u is above the lower limit u l , the second term in (21) is positive, thus b 1 = 1 and b 2 = 0 for the inequality to hold. If u is below the lower limit u l , then b 1 = 0 and b 2 = 1. The same holds analogously for the upper limit in (22), see Table 1.
With the discrete variables connected to the input ranges, it is now possible to give the output equation of the saturation as again with the slope s =

Composition
Systems can be connected in three different ways: series, parallel, and feedback. In general, any eMTI system can be connected to another eMTI system using one of the mentioned connection types. Whether the coupled system is within the class of eMTI systems depends on the subsystems and the connection type [26]. The following example shows a simple series connection where the overall system does not remain within the eMTI class.
Example 3.1. Consider two first order eMTI systems connected in series. The first system has one input u 1 and two outputs y 1 described bẏ The second system has a 2-dimensional input u 2 and one output y 2 described bẏ A series connection of the two subsystems can be represented as state space model Due to the term x 2 1 , it is no longer possible to model the overall system explicitly multilinear, but only by more general polynomial functions. See Appendix B for a detailed development and the parameter matrices F and G.
In contrast, if the systems are represented in iMTI form, the overall system will still be within this class, which is formally stated by the following proposition.  Proof. Given two model matrices H 1 and H 2 of the subsystems, the model matrix H of the overall system has a column dimension according to the total number of variables of the overall system. The row dimension takes the sum of the row dimensions of H 1 and H 2 plus the number of additional connecting equations, which is dependent on the connection type. The overall dimensions of H are given in subsection 2.3. The equations to append for different connection types are given next. 1. Series connection ( Figure 3) has the implicit equation 2. Parallel connection (Figure 4) has the implicit equations 3. Feedback connection ( Figure 5) has the implicit equations Example 3.2. If Example 3.1 is written in implicit form, it can be seen that the overall system remains in the iMTI class. The iMTI formulation of system 1 is and for system 2 Using (30) for series connection, the implicit form of the overall model is It is clear that the last two connection equations allow to reduce the model to See Appendix B for a more detailed development of the equations.
Equation (38) shows that in general, auxiliary variables are required in addition to the overall input and output variables and state variables of the subsystems. Symbolic methods for reduction of the number of auxiliary variables and thus, equations are known and provided by standard algorithms, e. g., as implemented for the Modelica language in OpenModelica. This paper does not focus on these, but on the structure of the iMTI models, which additional enable the use of tensor decomposition methods for reduction which will be discussed in the next section.

Tensors
Real-valued tensors X ∈ ℝ I 1 ×I 2 ×⋅⋅⋅×I n are multidimensional arrays well known from physics and various other application domains, e. g., in data sciences. Well established linear control engineering methods mostly rely on matrices, i. e., 2D tensors and also current extensions to the nonlinear domain, e. g., linear-parameter-varying models do this, too. But for multilinear models, matrices are less appropriate parameter formats than general tensors. This is because of the intrinsic tensor structure of multilinear monomials which can be seen best by the small example in Figure 6.
The contracted product ⟨⋅|⋅⟩ of a tensor F ∈ ℝ I 1 ×⋅⋅⋅×I n ×J 1 ⋅⋅⋅×J m and a tensor G ∈ ℝ I 1 ×⋅⋅⋅×I n is defined as a tensor The contracted product for tensors of the same size is called inner (or element-wise) product. This can be used to define a multilinear function with x ∈ ℝ n represented by the inner (elementwise) product of a parameter tensor F ∈ ℝ × n 2 , using the simple notation for tensor spaces as with the monomial tensor M(x) of the same dimension, which is the tensor analogon to (3), [13].

Tensor decomposition
The representation of the coefficients of a multilinear polynomial as tensor F as in (40) enables the application of tensor decomposition methods which have proven to be ex-tremely powerful in many other application domains. This paper focusses on the so called canonical polyadic (CP) decomposition, but all other decomposition methods are also possible in principle [14].
A tensor X can be decomposed in a sum of a minimum number of r outer products, where r is the so called tensor rank. All factors can be represented as matrices F i ∈ ℝ I i ×r for each dimension i of the original tensor F, which can be abbreviated as An element of X is then given by which can be seen as the generalization of the dyadic product for dimensions larger 2. Figure 7 shows how the red element x(3, 1, 2) of the tensor is computed as Determing the exact rank of certain tensors is known to be a relevant but hard problem in mathematics [7]. For engineering applications, the problem of finding a CP decomposition of a predefined rank as approximation of an original tensor F is more relevant. Minimizing the norm ⟨ E | E ⟩ of the error tensor still is a non-convex hard problem, but sub-optimal solutions can be computed by tools like Tensorlab [30] and used for approximation in applications like the one given in this paper in Section 6.

Explicit CP models
Any eMTI model can represented by contracted tensor products [26] with the parameter tensors F = ℝ × (n+m) 2×n and G = ℝ × (n+m )2×p . Whereas the parameter tensors can have arbitrary rank, monomial tensors are always rank 1, so for a state vector x ∈ ℝ n and an input vector u ∈ ℝ m the mono- mial tensor has the decomposition After CP decomposition of the parameter tensors as well as the monomial tensor, the eMTI model readṡ with the factor matrices F i , G i ∈ ℝ 2×r for all states and inputs i = 1, ⋅ ⋅ ⋅ , (n + m) and F ϕ ∈ ℝ n×r and G ϕ ∈ ℝ p×r . The next step is the most important one for applicability and large-scale models, because here the curse of dimensionality is broken by factorization. If a CP decomposition of the parameter tensor F with arbitrary rank r is available, the computation of the right hand side of the state space model (6)-(7) is possible by using the factors and the rank-1 structure of the monomial tensor, e. g., for the state equation bẏ The overall number of operations for each factor are 2r multiplications and n + m additions, together with n + m multiplications for the element-wise Hadamard product ⊛ and the final multiplication with the n × r matrix F ϕ these are (n + 2) r multiplications and r + n + m additions which is linear in all dimensions n and m of the variables whereas already the number 2 n+m of monomials has exponential complexity. Thus, full tensor computations are reasonable only for small dynamic order n and number m of inputs, whereas computations with decomposed representations are easily done up to very large dimensions. Similar arguments hold for the output equation (48).

Implicit CP models
An iMTI model (12) represented in tensor form reads and the same holds true for the representation of the inequality constraints As discussed, for larger systems only decomposed parameter tensors make sense, i. e., the general model is given by as mentioned in subsection 2.3, (50) can be included in (51) and therefore the decomposed version of (50) is not given. Computations with large scale models always have to be done on the factors like because even constructing the full parameter tensor would lead to enormous memory and runtime requirements. The factor matrix L ϕ ∈ ℝ N×r where the row dimension N is the number of inequality constraints and r is the decomposed tensor rank. All other factor matrices L i ∈ ℝ 2×r for i = 1, ⋅ ⋅ ⋅ 2n + m + p have a row dimension of 2. The next Section will show how this structure can be exploited further and how the storage needed for the model can nearly be reduced by another 50 %.

Normalization
In this section, a linear state and input transformation for eMTI models is given and a normalization of iMTI models in CP tensor representation are discussed.

Linear transformation
For linear state and input variable transformations of multilinear models (6) in continuous time given by the transformed state transition matrix is given by where diag denotes the diagonal matrix with elements a i and is a transformation matrix. Linear transformations allow a numerical preconditioning of models which is used in the application example. The linear transformation can as well be done in tensor representation.

Normalized implicit CP tensor models
The CP tensor H of the decomposed iMTI model from section 4 can be represented in a normalized form. Several norms, such as the euclidean norm (2-norm) or the manhattan norm (1-norm) could be used in principle, but the latter enables very efficient computations and thus is given priority in this paper [17]. In [13] it is proposed to normalize the columns of the factors of a CP tensor which is applied in [12] using the norm 1 condition.
A minimal representation of H by the factor H ϕ ∈ ℝ N×r and only vectors instead of matrices for the other factors thus is possible, because the other elements of the matrix H i can be computed by (58). Moreover, any CP decomposed tensor H of an iMTI model can be transformed in a minimal normalized representation with factor vectors and (2σ (H i (1, r)) − 1) .
Some remarks: -For computational efficiency, H ϕ is allowed to have arbitrary column vectors, -for all reasonable iMTI models, ||H i (:, k)|| 1 > 0 holds, -the Heaviside function (8) allows to construct a twovalued sign function, -all factor vectors h i can be stored in a matrix of dimension ℝ (2n+m+p)×r .
In the following, an example of CP tensor normalization is given by an example of an easy linear system relevant for the next section, too.

Example 5.1. A PI controller can be represented by the implicit equationsẋ
y − x − K P u = 0, with integral gain K I , proportional gain K P , input u, state x and output y. The CP iMTI model can be represented as a rank 5 CP tensor H with factor matrices for each term Because the third and the fourth columns of all h i are equal and H ϕ has only one parameter in the third row for the first state equation and one parameter in the forth row for the second state equation, it is obvious, that the rank can be reduced to 4 by storing the third and fourth columns as one: It is now possible to crosscheck, that by (53) the model (61) can be recovered To illustrate the normalization of the factor matrices from Example 5.1, the Figure 8

Application example
The previous section showed that the non-closedness of the eMTI is overcome by introducing the iMTI form. The following section presents a first showcase example about the usage of this new representation by connecting an electrolyzer with a buck converter powered by an infinite DC bus. These applications exemplify local energy systems, where several multiphysical dynamics are present. Three main aspects are demonstrated here. Firstly, the nonlinear model of the electrolyzer is multilinearized and compared with a typical linear approximation. Then this multiphysical model is connected with a buck converter, where a linear averaged model is assumed, since the focus in this paper is on the connection of the two models. Lastly, this paper derives a representation of the nonlinear saturation limits of the buck converter in an iMTI form. Thereby, the basic functionalities of the new implicit form are demonstrated and the prospective of large-scale multilinear energy systems modeling is shown.

Electrolyzer
The nonlinear model of the polymer exchange membrane (PEM) electrolyzer is based on a 46 kW model from literature [8]. Because in this example the connection to power electronics is crucial, also the electric dynamic behavior should be modeled [23]. Considering an equivalent electric circuit to model the electrolyzer with an RC element for each electrode, the ordinary differential equation (ODE) describing the activation overpotential v act on anode and cathode, respectively, is expressed as with the total electrolyzer current i ely in A [31]. The double layer capacity C DL , defined in (A.1), and activation resistance R act are often measured by a step function of the input current or voltage [10]. Since the base model of [8] does not include values for the capacity, simplifying assumptions will be made to include a value. The electric dynamic behavior is mainly caused by the double layer at the electrode-electrolyte interface. Due to the geometric arrangement of the electrode and electrolyte, the capacity can be estimated by a simple plate capacitor [29]. The total activation overpotential is mainly caused by the reaction at the anode (v act,an = v act ), wherefore the cathode is neglected and the resistance R act can be calculated as displayed in (A.2) [8]. Therefore, R act is already nonlinear, but in addition the temperature dependency of the electrode exchange current density j 0 is expressed by an exponential function [8]. The temperature influence on the open circuit voltage v ocv , defined by the Nernst-Equation, is also nonlinear (A.3). By adding up v ocv with the ohmic overpotential v ohm based on (A.4) and the activation overpotential, the total electrolyzer voltage can be described as Again, the temperature dependency of the protonic membrane conductivity σ mem is nonlinear and contains an exponential expression. For describing the thermal dynamics, an ODE is used as in with s 1 power (Ẇ), s 2 heat (Q) and s 3 enthalpy (Ḣ) streams, and the lumped thermal capacity C th [8]. The cooling heatQ cool is directly controlled by a PI controller, and described by a state space model aṡ with K P = −152.023, K I = −0.218 and the temperature difference as control deviation ΔT = T set − T. Because of the ODEs describing the thermal, electric and controller dynamics, the model has three states (n = 3). Since the electrolyzer is current-controlled, the input current i ely (m = 1) results into an output voltage v ely , which is set by the converter. Additionally, the generated hydrogen substance streamṅ H 2 in mol s −1 is considered as an output (p = 2), described by Farady's law iṅ containing the number of cells n c , Faraday's constant F and faradaic efficiency n F = 1. With the state and output vectors a nonlinear state space model can be derived.
Considering the broad operating range of electrolyzers, linearization of this nonlinear model at one operation point could be insufficient regarding the model accuracy. Thus, the projection multilinearization method from [16] was chosen, as it can be defined for the operating range of states, input and outputs. In that operating range the number of values per variable is N = 30 and the maximal multilinear order is m n = 2. The lower and upper bound for the multilinearization were estimated by the minimum and maximum values of the nonlinear model in the simulated range of 200 to 400 A, which corresponds to 50 to 100 % of the nominal load of the electrolyzer. This leads to the matrices F ∈ ℝ 3×16 and G ∈ ℝ 2×16 of the eMTI model. The eMTI model in the form of (6) and (7) can easily be transferred to implicit form, as described in (13) and (14) in subsection 2.3.
The eMTI model was decomposed by the MATLAB tensor toolbox, resulting in a reduced rank of r F = 9 and r G = 8. A linear transformation was applied to the model according to (54)-(56) for getting an operating region between 0 and 0.5. For numerical comparison, linearization at 300 A, corresponding to 75 % of the nominal load was done by the MATLAB function linmod.
At the operation point, the linear time-invariant (LTI) model fits the nonlinear dynamics better than the eMTI model, as it can be observed for v act in Figure 9 (c). After 10 min of electrolyzer operation at 75 %, the current jumps to nominal load, where the eMTI model fits better. The same behavior can be observed at the second jump to 200 A. Especially the states of the thermal dynamics T and x PI of the eMTI model show good accordance with the nonlinear model in Figure 9 (b) and (d). This example shows the benefits of the eMTI modeling approach for wider operating ranges compared to linearization.
Additionally, the computational time of the simulations are compared. The simulations are performed on a conventional laptop, while using the automatic solver and step-size selection within Simulink. The simulation time of the LTI model is the smallest (21 s), followed by the eMTI (26 s)  of the LTI model and the normalized decomposed tensors F, G of the eMTI model in double precision are used. The required memory of the eMTI model (888 bytes) is significantly higher than that of the LTI model (160 bytes). However, the advantages of decomposed tensors for MTI models will be noticeable for larger systems as the increase of elements is linear compared to an exponential increase for full tensors and for the LTI model. An example for this is given in Figure 10. The intersections of the graphs depend on the rank of the decomposed tensors. The memory requirements of the nonlinear model are not taken into account due to the significantly larger memory of the overall Simulink file. Therefore, a reasonable comparison is not possible here.

Buck converter
The buck converter steps down the voltage of 800 V provided by the infinite DC bus V DC . This paper takes an ideal model of a buck converter as depicted in Figure 11. At the load side of the converter, since the electrolyzer model provides the voltage v ely for a given current input i ely , a controlled voltage source is used to connect the electrolyzer to the DC bus. It is connected in series with a small resistance R S of 10 µΩ. The switching power electronic S, e. g., an IGBT, is modelled as an ideal switch. The diode, inductor L, and capacitor C, have no internal, or parasitic, re- sistances. The inductor is sized to 1 mH and the capacitor to 680 µF similarly to the filter design used in [9].
This paper assumes an average model of the switching, and thus, the buck converter excluding the controller can be described by the LTI state-space model where the duty cycle d is the input coming from the controller. This LTI model can be easily transformed into an iMTI model as outlined in (13)- (14) in subsection 2.3.

Local energy system
Finally, the coupled system of the buck converter and electrolyzer including a controller is analyzed. Figure 12 shows the closed control loop of the coupled system. The controller tracks the difference between the reference power P ref and the actual power P delivered to the electrolyzer. Using the tracking error, the duty cycle d is calculated with a PI controller, with gains K P = 5.0 ⋅ 10 −4 and K I = 1.0 ⋅ 10 −2 , and saturated to [0, 1]. The buck converter takes the duty cycle d and the output voltage of the electrolyzer v ely to provide the input current i ely to the electrolyzer.
As mentioned in section 3 the connection of eMTI systems can lead to an overall system, which is not in the eMTI class. Due to Lemma 3.1 the application example is formulated in implicit form, which guarantees the overall model to be within the iMTI class. The connection equations for Figure 11: Ideal buck converter scheme and connection to the electrolyzer.
the system are given as The iMTI representation of the controller can be written in CP tensor form where  (21) and (22). Special differential algebraic equations (DAE) solvers for descriptor models can be used for simulation of CP iMTI models like the one derived for the closed loop energy system here [4]. For this paper, the closed loop model is implemented in Simulink with CP eMTI submodels. Three different setups are compared, where the electrolyzer is either modeled as an LTI, eMTI or nonlinear system. The formulation of the electrolyzer w.r.t these three model classes is given in subsection 6.1. The rest of the system remains unchanged during the comparison shown in Figure 13.
A total simulation time of 30 minutes is used. During the simulation the reference power P ref is changed twice using a step function. The step is performed w. r. t. subsection 6.1. The initial reference power P ref until 10 minutes simulation time is related to 75 % of the nominal power of the electrolyzer. After the first step the reference power is increased by 1/3. After another 10 minutes PI (s) u Buck Converter  of simulation time the reference power is decreased by 2/3 of the initial power. Figure 13 gives an overview of the entire simulation regarding the time and operational range. Here, all models show comparable results.
A more detailed view of the results is given in Figure 14. An overshoot for the controller can be seen for all models after the first step. During the simulation from 10 to 20 minutes the eMTI model is outperforming the LTI model. Especially for an advanced simulation time after the first step, the deviation for the LTI model is large compared to the eMTI model.
The whole model can also be represented by a single tensor and inequalities where L tot appends the equalities to the inequalities given by L c . This can be achieved by expressing each equality by two opposing inequalities as given in (18). Although the number of inequalities gets larger, such a representation could lead to reductions because of the possible good low rank tensor approximations exploiting the internal structure of L tot .

Conclusion
This paper introduces iMTI models as a new class -which is interesting for large-scale hybrid complex engineering applications because of three reasons. The parameters can be represented as normalized tensors -which allow tremendous model reduction by modern decomposition methods. Nonlinear effects can be modeled -which is essential for the large signal behaviour of real world systems. Boolean dynamics is exactly modeled because its intrinsically multilinear -which is necessary to describe all discrete-valued subsystems, switches, and automata. Possible white-box modeling steps are shown by a small energy system example, including multilinearization, decomposition of the multilinear model into factors, normalization and composition back to overall models. Research on multilinear black-box and grey-box parameter identification is parallel ongoing. Future research activities will apply the iMTI form to modeling large power networks, where discrete-valued dynamics, such as the switches in power converters, become more present with the growth of renewables. This includes the quick composition of multiphysical energy systems models by using the implicit form. This will add a powerful tool to address the challenges arising in the European power system growing in complexity, for which efficient models are needed. with membrane thickness δ mem , standard protonic conductivity σ mem,std , and activation energy required for the proton transport in the membrane E pro [8].
For further description of the parameters and the thermal dynamics see [8].

Appendix B. Composition
This section provides the parameter matrices F and G of the eMTI model example in section 3. The state equation of system 1 is given bẏ