Implicit multilinear modeling

1.1 Motivation

Automation systems in general are hybrid, i. e., some signals are discrete-valued and some are continuous-valued [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.

Figure 1

Classes of functions.

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.

1.2 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.

1.3 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?

2.1 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.

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 B = { 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 ∈ 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

where x ¯ : = 1 − x ∈ R.

The following example of order n = 2, taken from [28], illustrates the idea.

2.2 Explicit multilinear models

With the monomial vector, an eMTI state space model can be written in continuous time

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.

Figure 2

Saturation function.

To model this as an eMTI model two discrete-valued variables are introduced

(9) b l = σ ( u − u l )

(10) b u = σ ( u − u u ) .

The output is then calculated by

(11) y = ( 1 − b l ) y l + b l ( 1 − b u ) ( ( u − u l ) s + y l ) + b u y u ,

with the slope s = y u − y l u u − u l of the middle part.

2.3 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 x ˙ ∈ R n , x ∈ R n , the inputs are u ∈ R m , the outputs are y ∈ R p and the model matrix is H ∈ R N × 2 2 n + 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.

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.

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.

4.1 Tensors

Real-valued tensors X ∈ R 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.

Figure 6

Monomial tensor of 3rd order.

The contracted product ⟨·|·⟩ of a tensor F ∈ R I 1 × ⋯ × I n × J 1 ⋯ × J m and a tensor G ∈ R I 1 × ⋯ × I n is defined as a tensor H = F G ∈ R J 1 × ⋯ × J m with elements

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 ∈ R n represented by the inner (elementwise) product of a parameter tensor F ∈ R × 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].

4.2 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 extremely 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 ∈ R 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 F 1 ( 3 , 1 ) F 2 ( 1 , 1 ) F 3 ( 2 , 1 ) + F 1 ( 3 , 2 ) F 2 ( 1 , 2 ) F 3 ( 2 , 2 ) + ⋯ + F 1 ( 3 , r ) F 2 ( 1 , r ) F 3 ( 2 , r ) by r · ( n − 1 ) multiplications and ( r − 1 ) summations.

Figure 7

CP Decomposition of a tensor with dimensions 5×3×4.

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.

4.3 Explicit CP models

Any eMTI model can represented by contracted tensor products [26]

After CP decomposition of the parameter tensors as well as the monomial tensor, the eMTI model reads

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 by

The overall number of operations for each factor are 2 r 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).

4.4 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 ϕ ∈ R 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 ∈ R 2 × r for i = 1 , ⋯ 2 n + 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 %.

5.1 Linear transformation

For linear state and input variable transformations of multilinear models (6) in continuous time 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.

5.2 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 ˜ ϕ ∈ R 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

Some remarks:

1. For computational efficiency, H ˜ ϕ is allowed to have arbitrary column vectors,

2. for all reasonable iMTI models, | | H i ( : , k ) | | 1 > 0 holds,

3. the Heaviside function (8) allows to construct a two-valued sign function,

4. all factor vectors h ˜ i can be stored in a matrix of dimension R ( 2 n + m + p ) × r .

To illustrate the normalization of the factor matrices from Example 5.1, the Figure 8 shows the original column vector of H u ( : , 3 ) with dashed line and the normalized column vector of H ˜ u ( : , 3 ) with solid line in blue. Red vectors belong to a factor H u ( : , k ) = 2 1 ′ to demonstrate normalization with arbitrary angle.

Figure 8

Normalization of a CP tensor factor to norm-1.

6.1 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 ( W ˙), s 2 heat ( Q ˙) and s 3 enthalpy ( H ˙) streams, and the lumped thermal capacity C th [8]. The cooling heat Q ˙ cool is directly controlled by a PI controller, and described by a state space model as

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 n ˙ H 2 in mol  s − 1 is considered as an output ( p = 2), described by Farady’s law in

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 ∈ R 3 × 16 and G ∈ R 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.

Figure 9

Electrolyzer model comparison for: (a) current input i ely , (b) temperature state T, (c) activation voltage state v act , and (d) cooling power state x PI .

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) and the nonlinear model (35 s). These results are expected as the LTI and eMTI are approximations of different complexity of the nonlinear model. In addition, the memory requirements are compared for the LTI and the eMTI models. For the comparison the matrices A, B, C, D 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.

Figure 10

Comparison of the memory consumption of an LTI model, eMTI model using full tensors and eMTI model using decomposed and normalized tensors with r F = 9 and r G = 8 as a function of number of states.

Figure 11

Ideal buck converter scheme and connection to the electrolyzer.

6.2 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, resistances. 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

Figure 12

Closed control loop.

6.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 the system are given as