Approach for the Mode Switching Problem in Piecewise Smooth Implicit

Multilinear IVPs

Torben Warnecke

1 a

and Gerwald Lichtenberg

2 b

Deutsches Elektronen-Synchrotron DESY, Notkestraße 85, 22607 Hamburg, Germany

Faculty of Life Sciences, Hamburg University of Applied Sciences, Ulmenliet 20, 21033 Hamburg, Germany

Keywords:

Hybrid Dynamical Systems, Multilinear Algebra, Differential-Algebraic Models, Descriptor Systems,

Switched Systems, Piecewise-Smooth Dynamical Systems.

Abstract:

This paper addresses the mode switching problem in piecewise smooth implicit multilinear initial value prob-

lems (IVPs), which are relevant for modeling hybrid dynamical systems like HVAC and power systems. Unlike

traditional switched systems with explicit mode descriptions, this work focuses on systems where mode in-

formation is implicitly encoded in binary-valued variables and switching conditions are deﬁned by inequality

constraints. The paper investigates the transversal motion discontinuities that occur when the system meets the

boundary surfaces deﬁned by these constraints. A method is presented to determine the discontinuous motion

by analyzing the total derivative of the inequality constraints. The modeling framework utilizes hybrid implicit

multilinear time-invariant (iMTI) functions and describes the system using inequality-constrained index-1

differential-algebraic equations (DAEs). The Jacobian matrices and thus the total derivatives can be estimated

algebraically, due to the use of multilinear functions. To handle the combinatorial complexity associated with

the binary variables during mode switching, the paper proposes using sparsity pattern analysis to identify and

solve sub-problems more efﬁciently. The presented method is applied to a two-point temperature-controlled

three-tank system, and simulations are performed using the MTI-Toolbox for MathWorks MATLAB.

1 INTRODUCTION

In recent researches implicit multilinear time invari-

ant (iMTI) models are of special interest since they

can describe fundamental dynamic relations in HVAC

and power systems, (Samaniego et al., 2024), (En-

gels et al., 2024), (Kaufmann et al., 2023), (Luxa

et al., 2022), (Pangalos et al., 2014). The multilinear

model class is a hyperclass of linear and binary mod-

els, while being a subclass of polynomial and general

nonlinear models, illustrated in ﬁgure 1 (Lichtenberg

et al., 2022).

The mode switching problem in piecewise smooth

implicit multilinear initial value problems (IVPs) is

important to simulate hybrid iMTI models reliably,

which gained attention due to the ability to model sat-

uration, switches and logical functions (Warnecke and

Lichtenberg, 2024a), (Lichtenberg et al., 2022). The

mode switching can be understood as the transver-

sal motion between the boundary surfaces of hy-

https://orcid.org/0009-0004-3037-8634

https://orcid.org/0000-0001-6032-0733

brid/switched iMTI models.

Usual switched systems have an explicit descrip-

tion of modes, switching conditions and transition

functions. They introduce multiple sets of equations

and constraints (Benzaouia, 2012) alongside multiple

sets of switching conditions either by equations, in-

equalities (Calvo et al., 2016) or subsets S ⊂ R

the signal space R

(Dieci and Lopez, 2011), as well

as transitions functions in between modes (Mehrmann

and Wunderlich, 2009).

Here a type of continuous-time hybrid iMTI mod-

els is of interest, described by implicit multilinear

DAEs, which only implicitly encode their mode in-

formation in binary-valued variables, while switching

conditions and transition functions are encoded in ad-

ditional inequality constraints. They include contin-

uous states x ∈ R

, input signals u ∈ R

, as well

as continuous- and binary-valued algebraic variables

y ∈ R

and z ∈ B

, with B = {true,false}.

The hybrid iMTI model framework shares struc-

tural similarities to the hybrid DAE model formula-

tion described in Appendix B of the Modelica Lan-

guage Speciﬁcation (Fritzson and Engelson, 1998).

Warnecke, T. and Lichtenberg, G.

Approach for the Mode Switching Problem in Piecewise Smooth Implicit Multilinear IVPs.

DOI: 10.5220/0013578200003970

In Proceedings of the 15th International Conference on Simulation and Modeling Methodologies, Technologies and Applications (SIMULTECH 2025), pages 349-356

ISBN: 978-989-758-759-7; ISSN: 2184-2841

349

linear

nonlinear

polynomial

binary

multilinear

Figure 1: Illustration of system classes.

Both approaches represent fully implicit formulations

of hybrid systems, lacking an explicit partitioning of

system dynamics for individual modes. However,

they differ in two key aspects: the iMTI framework

imposes a restriction to multilinear functions, and it

employs implicit inequality constraints in place of the

discrete equations used in Modelica.

Hybrid iMTI models offer two principal advan-

tages: First, they can reduce storage requirements

compared to general nonlinear switched systems.

This is due to the fact that the parameters and the

structure of both the equations and the inequality

constraints can be compactly represented through the

structure and parameter matrices S and φ, respec-

tively—eliminating the need for additional semantic

information, similar to the use of system matrices A,

B, C, and D in linear state-space models. Second,

unlike general nonlinear DAE formulations (such as

those described in the Modelica Language Speciﬁca-

tion), the constrained class of iMTI models allows for

algebraic derivation of Jacobian matrices with the di-

rect use of the structure and parameter matrices. This,

in turn, enables direct computation of total derivatives

of the model functions, as detailed in Chapter 2.1.

The publications solely focuses on the transver-

sal motion discontinuities of such IVPs and is limited

to index-1 DAEs without structural changes, e.g. a

change in the number of variables or variable-index

systems. Also sliding motions and their Filippov so-

lutions are not considered so far.

In (Calvo et al., 2016) an algorithm to determine

the discontinuous motion and solving the IVP for ex-

plicit piecewise smooth (switched) ODEs based of the

total derivative of the mode constraints is described.

Here, this idea is used to determine the implicitly de-

scribed discontinuity and solving the transversal mo-

tion problem of the piecewise smooth implicit multi-

linear IVPs. The method is described in chapter 3.

A block lower triangular (BLT) partitioning of the

incidence matrix is used to reduce complexity of the

mode switching problem, similar as in other tools, e.g.

Modelica (Sanz et al., 2014). The BLT partitioning

approach is described in chapter 4.

2 MODELLING FRAMEWORK

2.1 Hybrid Implicit Multilinear Models

The signal space

v(t) =







x(t)

u(t)

y(t)

z(t)







∈ R

× R

× B

(1)

of a hybrid implicit multilinear state-space model

consists of a continuous-valued state vector

x(t) ∈ R

, a continuous-valued input vector

u(t) ∈ R

, a continuous-valued algebraic signal

vector y(t) ∈ R

and a binary-valued algebraic

signal vector z(t) ∈ B

at continuous time t ∈ R.

Deﬁnition 2.1 (multilinear function). A multilinear

function is a function f

: R

→ R which is linear in

each individual variable separately, i.e. if all except

one arbitrary variable is held constant, it is a linear

function of this variable (Lang, 2002).

Deﬁnition 2.2 (hybrid multilinear function). A hybrid

multilinear function is a function f : H

→ R from a

hybrid space to a real number which gives the same

result as the corresponding multilinear function f

→ R if the binary variables are converted with the

standard injector

˜v











1 ∈ R if v

= true ∈ B

0 ∈ R if v

= f alse ∈ B

if v

∈ R.

(2)

SIMULTECH 2025 - 15th International Conference on Simulation and Modeling Methodologies, Technologies and Applications

350

Deﬁnition 2.3 (hybrid multilinear vector function).

A hybrid multilinear vector function is a function

f : H

→ R

which maps a hybrid space to a real

vector space where each element is given by a hybrid

multilinear function.

A continuous time hybrid implicit multilinear

time-invariant (hybrid iMTI) model is characterized

by a set of equations and inequality constraints that

represent an implicit system of inequality constrained

differential algebraic equations (DAEs)

f(v(t)) = 0, (3)

g(v(t)) ≤ 0, (4)

where f : R

× R

× B

→ R

n+p

and

g : R

× R

× B

→ R

are hybrid multi-

linear vector functions.

Deﬁnition 2.4 (Discontinuity surface of a hybrid

iMTI system). For the continuous-time hybrid dy-

namical system described by the equations (3) and

inequality constraints (4), the zeros of the element of

g(v(t)) span the discontinuity surfaces S. When a dis-

continuity surface S

is met at the event time t

, mean-

ing the function value of g

(v(t = t

)) = 0 is equal to

zero, the system can exhibit a discontinuous motion

and change its operation mode.

Deﬁnition 2.5 (Operation mode of a hybrid iMTI sys-

tem). Different from classical switched system ap-

proaches the operation mode of the system is not

explicitly stated nor governed by individual sets of

equations and constraints. Instead, the operation

mode is implicitly described by the combination of

the binary-algebraic variables z. The Domain of each

mode is constrained by the inequality constraints (4).

Moreover, not all combinations of z must be accessi-

ble and therefore the number n

mode

of possible opera-

tion modes is less or equal to the number of possible

combinations of z ∈ B

mode

≤ 2

. (5)

The multilinear functions of a hybrid iMTI model

can be written in a normalized canonical polyadic de-

composition (J

ores et al., 2022):

(v(t)) =

∑

k=1

2n+m+p+q

∏

i=1

(1− | S

| +S

(t)),

(6)

(v(t)) =

∑

k=1

2n+m+p+q

∏

i=1

(1− | S

| +S

(t)),

(7)

where all model information is stored in the structure

matrix S ∈ R

(2n+m+p+q)×R

and parameter matrices

∈ R

(n+p)×R

and Φ

∈ R

q×R

Deﬁnition 2.6 (Jacobian matrix of a hybrid multilin-

ear vector function). The Jacobian matrix Jf(v(t))

is a matrix arrangement of the ﬁrst derivatives of

a vector function f(v(t)) ∈ R

(n+p)

with the vari-

ables v(t) ∈ H

, where an element J f

(v(t)) rep-

resents the ﬁrst derivative

∂ f

(v(t))

∂v

of the element f

with respect to the variable v

. The Jacobian ma-

trix Jf ∈ R

(n+p)×(2n+m+p+q)

of the hybrid multilinear

vector function f can be directly computed using the

normalized CP decomposed representation, as shown

in (Kaufmann et al., 2023). Each element

J f

(v(t)) =

∑

k=1

∏

i ∈ N

(1− | S

| +S

(t))

(8)

can be computed at time t with the signal vector

v(t), the structure matrix S and the parameter matrix

∈ R

(n+p)×R

, where N = {1...(2n+m+ p +q)}\b

is the set of indices of the variables excluding the in-

dex b of the variable v

(t). Similar for g(v(t)) ∈ R

Deﬁnition 2.7 (Total derivative of a hybrid multilin-

ear vector function). The total derivative

Df(v(t)) = Jf(v(t))

v(t)

∈ R

n+p

(9)

of a hybrid multilinear vector function

f(v(t)) ∈ R

(n+p)

with v(t) ∈ H

(2n+m+p+q)

time t can be computed from its Jacobian matrix

Jf ∈ R

(n+p)×(2n+m+p+q)

at time t and the time-

derivative

v(t)

∈ H

(2n+m+p+q)

of the signal vector

v(t) (using the standard injector). The total deriva-

tive describes the evolution of the function values

over time.

Remark 2.1. The time derivative of the binary vari-

ables is assumed to be zero, since a switch in its value

can be understood as switching the systems equations.

So they are treated as constant parameters of the dif-

ferent functions. Or differently one could state, that

the time derivatives of z directly before and after the

switch is zero.

2.2 Example: 2-Point

Temperature-Controlled Three

Tank System

The system is illustrated in ﬁgure 2. It consists of 3

water tanks sharing one inlet ﬂow of temperature u

The water tanks transfer heat between each other, as

well as losing heat to the environment with the am-

bient temperature u

. If multiple valves (z

, z

and

) are opened at the same time the inlet ﬂow will be

shared uniformly between those. The three 2-Point

Approach for the Mode Switching Problem in Piecewise Smooth Implicit Multilinear IVPs

351

controllers are designed to maintain the temperatures

of the tanks (x

, x

and x

) between lower (u

, u

and

) and upper bounds (u

, u

and u

Figure 2: Schematic of the hybrid temperature-controlled

three tank system.

The system can be modelled by the following set

of implicit inequality-constrained multilinear DAEs:

(cm ˙x

+ γx

− α

+ x

) − α

+ c f

− u

= 0,

(10)

(cm ˙x

+ γx

− α

+ x

) − α

+ c f

− u

= 0,

(11)

(cm ˙x

+ γx

− α

+ x

) − α

+ c f

− u

= 0,

(12)

− z

+ z

− 1 = 0, (13)

− x

)(1 − z

+ z

) ≤ 0, (14)

− x

)(1 − z

+ z

) ≤ 0, (15)

− x

)(1 − z

+ z

) ≤ 0, (16)

− u

)(1 − z

+ z

) ≤ 0, (17)

− u

)(1 − z

+ z

) ≤ 0, (18)

− u

)(1 − z

+ z

) ≤ 0, (19)

where γ = (α

+ 2α

), and the parameters are

the heat capacity c = 4200kJ/kgK and mass m =

5kg of water in each tank, the ambient heat

transfer α

= k

= 125W/K and heat transfer

= k

= 125W /K in between the tanks, and

valve coefﬁcient f

= 0.09kg/s. The model is

feasible for



(t) u

(t)



⊺

≤ x(t) and

x(t) ≤



(t) u

(t)



⊺

3 METHOD

When one or more discontinuity surfaces S are met,

meaning the function values of one or more inequality

constraints g

(v(t = t

)) = 0 are equal to zero, where

is the event time and g

are the according elements

of the vector function g, one needs to determine how

the system will behave.

Usually a hybrid implicit multilinear system is ca-

pable of different motions:

• staying in the previous mode, so the elements of

will continuously go back to negative values

and the binary algebraic variables won’t change

• sliding along a boundary, where some elements of

will stay at zero (Filippov solutions, (Piiroinen

and Kuznetsov, 2008))

• change its mode (transversal motion), where the

binary variables will change, so that either the

function values g

discontinuously switch to neg-

ative values or continuously go back to negative

values after the switch

• a combination of transition and sliding.

When a boundary is met, which can be nu-

merically detected by event detection algorithms

implemented in ODE solvers, one needs to solve the

inequality constrained system of DAEs at the event

time t

, with testing possible combinations of z until

a solution for the discontinuous motion problem of

the dynamical system is found.

Solution of the Inequality Constrained System of

Equations at the Event Time t

The unknown state derivatives

x(t = t

) and

continuous-valued algebraic variables y(t = t

) can be

computed by solving the inequality constrained sys-

tem of equations (3) at the event time t

for differ-

ent combinations of z, where the solution must sat-

isfy (4). A solution can be efﬁciently obtained by

the method described in (Warnecke and Lichtenberg,

2024b) where the problem is decomposed in sub-

problems and explicit solutions.

It is important to mention, that the solution to

the inequality constrained system of equations at the

event time t

must not be a solution for the dynamical

problem of the inequality constraint system of DAEs.

Usually at the time t

there is more than one solution

to the inequality constrained system of equations,

since the current system state is already one of them.

Solution of the Discontinuous Motion Problem.

Similar to the solutions of the inequality constrained

system of equations, the solutions for the dynami-

cal problem of the inequality constrained system of

DAEs must satisfy (3) and (4). If all elements of the

vector function

g(v(t = t

)) < 0 (20)

are negative the solution is also a solution for the

dynamical system, and it will undergo a transver-

sal motion, where the according variables z repre-

sent the future mode of the system. If some values

of g(v(t = t

)) are equal to zero:

(v(t = t

)) = 0 (21)

SIMULTECH 2025 - 15th International Conference on Simulation and Modeling Methodologies, Technologies and Applications

352

with l as the indices for the according elements of

g(v(t = t

)), the according elements Dg

(v(t = t

)) of

the total derivative Dg(v(t = t

)) of g(v(t = t

)) must

be equal to or smaller than zero:

(v(t = t

)) ≤ 0, (22)

so that the solution is also a solution for the dynami-

cal system. If Dg

(v(t = t

)) = 0 the system will slide

along the boundary. If the combination of the future

and current binary variables z are not similar the sys-

tem undergoes a transversal motion.

Remark 3.1. If for all possible combinations of z the

function values g(v(t = t

)) ≤ 0 are smaller or equal

to zero and the total derivatives Dg

(v(t = t

)) > 0

of boundaries g

(v(t = t

)) = 0 is bigger than zero,

the behavior can be described as zeno behavior, since

even when the mode changes, the system will immedi-

ately be pushed towards the boundary.

To compute Dg(v(t = t

)) the missing time-

derivatives

and

of the signal vector v(t) can be

computed with the total derivative Df(v(t)) of the vec-

tor function f(v(t)) of the equations (3), with the as-

sumption that the time-derivative of the binary-valued

algebraic signals is equal to zero right after the switch.

Since the value of the vector function is constant and

equal to zero at every time t, the total derivative

Df(v(t)) = Jf(v(t))







x(t)

u(t)

y(t)

z(t)







= Jf(v(t))







x(t)

u(t)

y(t)







= 0

(23)

is also equal to zero at every time t as well.

(23) can be rearranged for calculating

x(t) and

y(t) as the solution of the linear system of equations

x,y

(v(t))



x(t)

y(t)



= −Jf

x,u

(v(t))



x(t)

u(t)



, (24)

if Jf

x,y

(v(t)) is invertible, where Jf

x,u

(v(t)) repre-

sents the columns of the Jacobian matrix according

to the states x(t) and inputs u(t) and Jf

x,y

(v(t)) rep-

resents the columns of the Jacobian matrix according

to the state derivatives

x(t) and continuous-valued al-

gebraic signals y(t) at time t. The derivatives

u(t)

can be computed from the given input signal u(t), e.g.

when using linear interpolation for a given set of input

points as the slopes of the linear piece-wise function.

With

x(t),

y(t) and

u(t) the total derivative

Dg(v(t)) = Jg(v(t))







x(t)

u(t)

y(t)







(25)

of the inequality constraints g(v(t)) can be computed.

Remark 3.2. If one is only interested in one solution,

and it can be any of the possible solution, it would be

easy to iterate through all possible combinations until

an appropriate solution is found.

Remark 3.3. Numerical precision can be crucial for

solving the discontinuous motion problem, because

the possible solutions differ drastically around the

function values of zero and slight inaccuracy can

change the result, especially regarding the inequality

constraints and binary variables. To solve such nu-

merical issues normalization techniques for the func-

tion values could be used, e.g. by normalizing both,

the signal and parameter space. Also, one could think

about unilateral tolerances of the numerical proce-

dure.

Still the problem remains of high combinatorial

complexity with 2

possible combinations of z. For

most systems it would be reasonable to ﬁrst try the

solutions with the smallest hamming distance. Also,

in most IVPs the problem can be reduced by sparsity

pattern analysis, as described in the following chapter.

4 PROBLEM REDUCTION BY

SPARSITY PATTERN ANALYSIS

The sparsity pattern of a system of equations and

constraints can be represented by the incidence

matrix P of the IVP. By exploiting sparsity, e.g. by

reordering of the sparsity pattern, one can reduce the

complexity of the discontinuous motion problem. A

sparsity pattern P of a (hybrid) iMTI model can be

calculated as described in (Warnecke and Lichten-

berg, 2024b). By ordering the sparsity pattern of the

unknown variables in a block lower triangular matrix

P distinct sub-sets and sub-problems of the system

of inequality constrained DAEs can be identiﬁed

(Baharev et al., 2019). Different heuristics can be

used for the NP-complete problem (Vassilevska

and Pinar, 2004), e.g. by using algorithms like the

Dulmage–Mendelsohn decomposition (Dulmage and

Mendelsohn, 1958), as in (Pothen and Fan, 1990), or

simple heuristics as in (Pinar and Heath, 1999).

Deﬁnition 4.1 (sub-set). Sub-sets are sets of equa-

tions and constraints of a system of implicit equa-

tions and constraints, which do not have any simi-

lar unknown variables (unknown regarding the IVP)

and therefore can be solved independent of each other

(Warnecke and Lichtenberg, 2024b).

The index set A denotes the set of indexes ac-

cording to the equations and constraints of the sub-set

Approach for the Mode Switching Problem in Piecewise Smooth Implicit Multilinear IVPs

353

(rows of the sparsity pattern

P). The index set B de-

notes the set of indexes according to the variables of

the sub-set (columns of the sparsity pattern

P). In a

block lower triangular sorted sparsity pattern

P, sub-

sets are diagonal blocks which do not have any over-

lapping blocks to other diagonal blocks. The elements

ˆp

of the rows and columns adjacent to such a block

must only consist of zeros:

ˆp

= 0, (26)

for i /∈ A and k ∈ B, as well as i ∈ A and k /∈ B.

Deﬁnition 4.2 (sub-problem). A sub-problem is a set

of equations and constraints, which is part of a sub-

set of a system of implicit equations and constraints

(Warnecke and Lichtenberg, 2024b).

The index-set C denotes the set of indexes ac-

cording to the equations and constraints of the sub-

problem (corresponding to rows in the sparsity pat-

tern

P). The index-set D denotes the set of indexes

according to the variables of the sub-problem (corre-

sponding to columns of the sparsity pattern

P). Un-

like sub-sets, different sub-problems can have similar

(unknown) variables.

In a block lower triangular sparsity pattern

P, sub-

problems are diagonal blocks which can have (lower

and left) overlapping blocks to other diagonal blocks,

but the elements ˆp

of the right and upper parts of the

rows and columns adjacent to the sub-problem must

only consist of zeros:

ˆp

= 0, (27)

for i < c, where c ∈ C and ∀

x ∈ C

c ≤ x, and k ∈ D, as

well as i ∈ C and k > D, where d ∈ D and ∀

x∈D

d ≥ x.

A sub-problem determines the unknown variables

according to its columns of the sparsity pattern, by

solving the set of implicit equations and constraints

according to its rows, but can be dependent on previ-

ously solved variables of other sub-problems.

For the example model the sorted sparsity pattern

P =







1 1 0 0 0 0 0 0 0 0

0 0 1 1 0 0 0 0 0 0

0 0 0 0 1 1 0 0 0 0

0 1 0 1 0 1 1 0 0 0

0 1 0 0 0 0 1 1 0 0

0 0 0 1 0 0 1 0 1 0

0 0 0 0 0 1 1 0 0 1







, (28)

regarding the unknown variables of the IVP (

x, y and

z) reveals one overall sub-set with 7 sub-problems.

The sub-problems can be computed in series, where

instead of 2

= 64 overall combinations of binary

variables, the number of possible combinations that

need to be tested reduces to 3 × 2

= 12, since they

are independent of each other.

With the use of the sub-problems and sub-sets, the

combinatorial complexity of the discontinuous mo-

tion problem can be reduced. In ﬁgure 3 the reduced

brute force method is shown.

denote a combination of the binary variables z,

where z

i+1

denote the next possible combination (e.g.

with the smallest hamming distance to z

). The values

of the binary variables z

which belong to a previous

sub-problem in each sub-set of the events g

are ﬁxed,

since they cannot be inﬂuenced by the events, thus

reducing the possible combinations of z.

5 SIMULATION OF THE

EXAMPLE

The three tank system is modelled and simulated

using the MTI-Toolbox (Lichtenberg et al., 2024)

for MATLAB, where the method is incorparated.

The smooth sections of the IVP are calculated us-

ing the ode15i-solver with a BDF integration scheme

(Shampine, 2002). If a boundary surface is met, lo-

cated by the event detection algorithm of the solver,

integration stops and the procedure from ﬁgure 3 is

executed. When a solution for the transversal motion

has been found, integration continues with the new

values of z.

The example is simulated for 1000 sec-

onds, with the intial states x(t = 0) =



37 30 40



⊺

and constant input signals

u(t) =



0 80 30 40 25 35 35 45



⊺

∀t. In

ﬁgure 4 the trends of the states x, valve positions z

and z

and algebraic signal y

can be seen. The

simulation has taken 1.19 seconds. It can be seen

that the system has state-continuity, which means

that the states cannot jump in mode transitions and

are smooth in between, while the algebraic variables

y can be discontinuous (jump) in between mode

transitions. Even if the continuous-valued algebraic

variables y ⊆ R are in general restricted to a subset

Y of the real number space R, in this example the

variable y

is restricted to a subset of the natural

number space {1,2,3} ⊂ N, since the signal y

exclusively dependent on z

, z

and z

In the phase diagram in ﬁgure 5 it can be seen, that

the trend of the states has corners at the discontinuity

surfaces marked by the black squares. More over the

diagram reveals a period-2 limit cycle.

SIMULTECH 2025 - 15th International Conference on Simulation and Modeling Methodologies, Technologies and Applications

354

Figure 3: Brute force procedure for solving the discontinuous motion problem when sub-sets and sub-problems are known.

Figure 4: Time signals of the simulations results of the three

tank model.

6 CONCLUSION

The publication focuses on the transversal motion dis-

continuities of piecewise smooth implicit multilinear

initial value problems. Unlike usual switched systems

with explicit descriptions of modes, switching condi-

tions, and transition functions, the hybrid iMTI mod-

els implicitly encode mode information in binary-

valued variables, with switching conditions and tran-

Figure 5: Phase diagram of the simulations results of the

three tank model.

sitions deﬁned by inequality constraints.

A method has been introduced to determine

the implicitly described discontinuity and solve the

transversal motion problem by adapting an idea from

explicit piecewise smooth ODEs based on the total

derivative of mode constraints.

A 2-point temperature-controlled three-tank sys-

tem is used to illustrate the application of the hybrid

iMTI modeling framework, as well as the method to

solve the piecewise smooth IVP. The method for han-

dling mode switching when discontinuity surfaces are

met involves solving an inequality-constrained sys-

Approach for the Mode Switching Problem in Piecewise Smooth Implicit Multilinear IVPs

355

tem of equations at the event time for different com-

binations of binary variables and checking conditions

based on the function values and their total derivatives

to determine the system’s behavior.

To address the combinatorial complexity arising

from the binary variables, it is proposed to use spar-

sity pattern analysis to identify sub-sets and sub-

problems within the system, which can be solved

more efﬁciently. The simulation results of the exam-

ple show the applicability of the method. The sys-

tem’s behavior at discontinuity surfaces can be seen

and reveal a limit cycle.

REFERENCES

Baharev, A., Neumaier, A., and Schichl, H. (2019). A

manifold-based approach to sparse global constraint

satisfaction problems. Journal of Global Optimiza-

tion, 75(4):949–971.

Benzaouia, A. (2012). Introduction to switched systems. In

Saturated Switching Systems, pages 77–94. Springer.

Calvo, M., Montijano, J., and R

andez, L. (2016). Algo-

rithm 968: Disode45: A matlab runge-kutta solver for

piecewise smooth ivps of ﬁlippov type. ACM Trans-

actions on Mathematical Software, 43:1–14.

Dieci, L. and Lopez, L. (2011). Fundamental matrix solu-

tions of piecewise smooth differential systems. Math.

Comput. Simul., 81:932–953.

Dulmage, A. L. and Mendelsohn, N. S. (1958). Coverings

of bipartite graphs. Canadian Journal of Mathematics,

10:517–534.

Engels, M., Lichtenberg, G., and Knorn, S. (2024). An

approach to parameter identiﬁcation for boolean-

structured multilinear time-invariant models*. In 2024

European Control Conference (ECC), pages 2077–

2082.

Fritzson, P. and Engelson, V. (1998). Modelica—a uni-

ﬁed object-oriented language for system modeling

and simulation. In ECOOP’98—object-oriented pro-

gramming: 12th European conference Brussels, Bel-

gium, July 20–24, 1998 proceedings 12, pages 67–90.

Springer.

ores, N., Kaufmann, C., Schnelle, L., Y

nez, C., Pangalos,

G., and Lichtenberg, G. (2022). Reduced cp represen-

tation of multilinear models. Proceedings of the 12th

International Conference on Simulation and Modeling

Methodologies, Technologies and Applications.

Kaufmann, C., Garcia, D. C., Lichtenberg, G., Panga-

los, G., and Y

nez, C. C. (2023). Efﬁcient lin-

earization of explicit multilinear systems using nor-

malized decomposed tensors. IFAC-PapersOnLine,

56(2):7312–7317.

Lang, S. (2002). Algebra. Graduate Texts in Mathematics,

volume 211, chapter Matrices and Linear Maps §4.

Determinants, page 511. Springer New York.

Lichtenberg, G. et al. (2022). Implicit multilinear modeling.

at - Automatisierungstechnik, 70(1):13–30.

Lichtenberg, G., Uhlenberg, E., Warnecke, T., Cateriano,

C., Engels, M., Samaniego, L., Schnelle, L., and

Kaufmann, C. (2024). MTI-Toolbox. http://mti.

systems.

Luxa, A. et al. (2022). Multilinear modeling and simula-

tion of a multi-stack PEM electrolyzer with degrada-

tion for control concept comparison. Proceedings of

the 12th International Conference on Simulation and

Modeling Methodologies, Technologies and Applica-

tions.

Mehrmann, V. and Wunderlich, L. (2009). Hybrid sys-

tems of differential-algebraic equations - analysis and

numerical solution. Journal of Process Control,

19(8):1218–1228. Special Section on Hybrid Sys-

tems: Modeling, Simulation and Optimization.

Pangalos, G., Eichler, A., and Lichtenberg, G. (2014). Hy-

brid multilinear modeling and applications. Advances

in Intelligent Systems and Computing, page 71–85.

Piiroinen, P. and Kuznetsov, Y. (2008). An event-driven

method to simulate ﬁlippov systems with accurate

computing of sliding motions. ACM Transactions on

Mathematical Software, 34.

Pinar, A. and Heath, M. T. (1999). Improving performance

of sparse matrix-vector multiplication. Proceedings of

the 1999 ACM/IEEE conference on Supercomputing.

Pothen, A. and Fan, C.-J. (1990). Computing the block tri-

angular form of a sparse matrix. ACM Trans. Math.

Softw., 16(4):303–324.

Samaniego, L., Uhlenberg, E., Kaufmann, C., Engels, M.,

Pangalos, G., Cateriano Y

nez, C., and Lichtenberg,

G. (2024). An approach to multi-energy network mod-

eling by multilinear models. pages 1522–1527.

Sanz, V., Urquia, A., and Casella, F. (2014). Improving

efﬁciency of hybrid system simulation in Modelica.

volume 2014.

Shampine, L. (2002). Solving 0=f(t,y(t),y ’ (t)) in matlab.

Journal of Numerical Mathematics, 10.

Vassilevska, V. and Pinar, A. (2004). Finding nonoverlap-

ping dense blocks of a sparse matrix. Lawrence Berke-

ley National Laboratory.

Warnecke, T. and Lichtenberg, G. (2024a). Hybrid implicit

multilinear modeling of complex hvac systems. In

Wagner, G., Werner, F., and De Rango, F., editors,

Simulation and Modeling Methodologies, Technolo-

gies and Applications, pages 79–104, Cham. Springer

Nature Switzerland.

Warnecke, T. and Lichtenberg, G. (2024b). Hybrid im-

plicit multilinear simulation using difference algebraic

equations reordering by sparsity patterns*. In Interna-

tional Conference on Control, Decision and Informa-

tion Technologies 2024, 2024 10th International Con-

ference on Control, Decision and Information Tech-

nologies (CoDIT), pages 2578–2583.

SIMULTECH 2025 - 15th International Conference on Simulation and Modeling Methodologies, Technologies and Applications

356