A distributional solution framework for linear hyperbolic PDEs coupled to switched DAEs

A distributional solution framework is developed for systems consisting of linear hyperbolic partial differential equations and switched differential-algebraic equations (DAEs) which are coupled via boundary conditions. The unique solvability is then characterize in terms of a switched delay DAE. The theory is illustrated with an example of electric power lines modelled by the telegraph equations which are coupled via a switching transformer where simulations confirm the predicted impulsive solutions.


Introduction
In this paper, we develop a rigorous solution theory for systems where a linear hyperbolic partial differential equation (PDE) is coupled with a switched differentialalgebraic equation (DAE) via boundary conditions (BC), see Fig. 1 as an overview.
Such systems occur, for example, when modelling power grids using the telegraph equation [8] including switches (e.g. induced by disconnecting lines), water flow  [13,14], supply chain models including processor breakdown [1,7], district heating systems with rapid consumption changes [5] and blood flow with simplified valve models in the heart [15]. Similar to [3,11] the closed-loop setting illustrated in Fig. 1 can include general network structures. In this coupled system, the output of the switched DAE provides the boundary condition for the PDE and the boundary values of the PDE serve as input to the DAE. Solutions of switched DAEs in general contain jumps and derivatives thereof, e.g. derivatives of Dirac impulses [17,19]; hence, the solution concept of the PDE has to be extended to allow for jumps, Dirac impulses and their derivatives at the boundary. In particular, this is a wider class compared to the solutions of small bounded variation, e.g. used in [4] where a nonlinear hyperbolic PDE is coupled to an ODE. Similarly, in [2,12], the investigations of switched linear PDEs with source terms are restricted to solutions with bounded variation. In [16], Dirac impulses are introduced at the position of an interface of nonlinear PDEs. A more general appearance of Dirac impulses is allowed in [6,23] for a partially linear system. A detailed overview of existing solution concepts for linear hyperbolic PDEs is given in [10].
Since arbitrary high derivatives of Dirac impulses can occur as solutions of switched DAEs, the aforementioned approaches are not suitable to handle the coupled systems studied here. Indeed, our first main contribution is a suitable extension of the 1D piecewise-smooth distributional solution framework (developed to handle switched DAEs in [17,18]) to a 2D piecewise-smooth distributional solution framework. This new class of solutions is large enough to allow distributional solutions such as Dirac impulses and their derivatives. At the same time, its elements are regular enough to allow a trace evaluation on the boundaries of the domain.
Towards our main existence and uniqueness result for solutions of the coupled system, we also establish a relationship between the solutions of the coupled systems and the solution of a switched delay DAE. For the latter, we generalize a recent existence and uniqueness result for delay DAEs in [20]. This paper is structured as follows.
After a detailed description of the coupled system (including an example of a simple power network), we review the classical solution theory of linear hyperbolic PDEs in Sect. 3 and the solution theory of switched DAEs in Sect. 4. The novel distributional solution framework for linear hyperbolic PDEs is introduced in Sect. 5 which is then used in Sect. 6 to establish a link between the coupled system and the solutions of a switched delay DAE (Theorem 23). Finally, we establish an existence and uniqueness result for general switched delay DAE (Theorem 24) and can conclude our main result about the existence and uniqueness of solutions of the coupled system (Corollary 26). We illustrate the results by numerical simulations of the power grid example.

System class
We consider linear hyperbolic PDEs on a bounded interval of the form: where a, b, t 0 ∈ R with a < b, u : [t 0 , ∞)×[a, b] → R n , n ∈ N, is the n-dimensional vector of unknowns of the PDE, A ∈ R n×n and y P : [t 0 , ∞) → R ν , ν ∈ N is the νdimensional output of the PDE depending on u ab (t) := (u(t, a) , u(t, b) ) ∈ R 2n and C P ∈ R ν×2n . The boundary conditions (BC) of the PDE have the form: where P ∈ R n×2n and y D : [t 0 , ∞) → R n is the output of the switched DAE with the m-dimensional vector of unknowns w : [t 0 , ∞) → R m , m ∈ N, the switching signal σ : R → {1, 2, . . . , N }, N ∈ N, and E ξ , H ξ ∈ R m×m , B ξ ∈ R m×ν , f ξ : [t 0 , ∞) → R m , C D ξ ∈ R n×m for each ξ ∈ {1, 2, . . . , N }.
The coupled system (2.1) has to be equipped with initial conditions where u t 0 : [a, b] → R n and w t 0 ∈ R m .

Remark 1
We would like to stress that the coupling structure in (2.1) is quite general; in fact, arbitrary networks whose edges represent PDEs and whose nodes represent (switched) DAEs which couple the different PDEs are covered. Consider, for example, a network as illustrated in Fig. 2a, where on each edge E the quantity u E is governed by a linear PDE u E t + Au E x = 0. At each node s, algebraic and/or differential conditions combine possible internal states w s with certain boundary values q s of the connected u E , i.e. E s σẇ s = H s σ w s + B s σ q s + f s σ . This setup can be rewritten in the form (2.1) by first rescaling the spatial domain (which simply modifies the matrices A E by a constant multiple) so that all PDEs are defined on the same interval and can be viewed as a single PDE where the new Fig. 2 An example of a network consisting of four nodes and five connecting edges and how to reduce it to one node and one edge (loop) network which still has all the features of the original network unknown u consists of the unknowns u E of each edge stacked over each other. (The A-matrix then is a block diagonal matrix.) In a similar way, the switched differentialalgebraic equations for each node can be stacked over each other (resulting in block diagonal coefficient matrices) and can then be viewed as single switched DAE, see Fig. 2b. A similar reduction is used in [3,11]. This method is also used in the specific example of a simple power grid discussed in Sect. 2.3.

Ill-posed coupling
As mentioned in the introduction, coupling a PDE to a switched DAE is challenging due to the need to consider distributional boundary terms. However, there is another challenge due to the additional algebraic constraints present in a DAE. The following example illustrates that even if the PDE is well-posed (i.e. it has a unique solutions for all initial and boundary conditions) and the DAE is well-posed (i.e. it has unique solutions for all consistent initial values and for all sufficiently smooth input signals), the coupled system is ill-posed.

Example 2 Consider the scalar advection equation given by
and output y P (t) := v(t, 0). It is well known that there exists a unique solution to (2.3) for any sufficiently smooth v 0 (·) and b(·). This PDE will be coupled to the following (nonswitched) DAE given by with input q(·) and output y D (t) := [ 0 1] w z . It is easily seen that for all q(·) and all initial values for w(0) there exists a unique solution. Now let (2.3) be coupled with (2.4) via b(t) = y D (t) and q(t) = y P (t). Since by construction b(t) = y P (t), this coupling leads to q(t) = y D (t) = z(t), which means that the DAE changes to However, in this new DAE, the variable z is completely free, which also means that the boundary of the PDE is completely free, i.e. the coupled system does not have a unique solution anymore. For this simple example, it is easy to see why well-posedness of the coupled system is lost: The PDE has as an "output" the same value as the provided boundary value; hence, both values are not independent and cannot be coupled via a DAE in an arbitrary way. In fact, just choosing a different output for the DAE (e.g. y D = [ 1 0] w z results in a well-posed coupled system.
In general, it will not be so apparent as in the above example which couplings lead to a well-posed system and it is therefore necessary to better understand the underlying solution structure.

Power grid example
Consider the simple electrical power grid illustrated in Fig. 3.
Each line is modelled by the telegraph equation of the form: where t ≥ 0, x ∈ [0, ], I (t, x) stands for the current, V (t, x) is the voltage and the constants L and C are inductance and capacitance, respectively. In particular, each line k has a position-dependent current I k and voltage V k . By appropriate scaling of the coefficients in the telegraph equation, we can assume that all PDEs are defined on the common domain T × X = [0, ∞) × [a, b]. At the nodes, there is a coupling between corresponding boundary values, where the "outputs" of the telegraph equations are the boundary currents I k for each line k. A generator is located at the first node, where we assume an externally given voltage. This algebraic constraint is modelled by the algebraic relations: where v G (·) is the externally given (time-varying) voltage of the generator together with the boundary condition V 1 (·, a + ) = y G D . On the consumption nodes, all voltages are assumed to be equal and we assume that the consumption is modelled as a simple Ohm's resistance, i.e. the sum of the (directed) currents at the boundary of the lines is proportional to the voltage at the node, and this is modelled by the DAEs where R 24 , R 34 > 0 are the resistive loads. Further, we impose the boundary conditions Finally, the switching transformer node is governed by the electric circuit given in Fig. 4.
The switch-independent equations governing this switching transformer node are as follows: where κ 12 , κ 13 > 0 are the voltage amplification constants. Note that, in this example, we use amplifiers only for the voltage values, so the power grid example is a simplified model. If the switch connects line 1 and 2, then the following two algebraic constraints hold together with the output y 1 D = v 12 and, otherwise, , then the rules governing the switching transformer node can be compactly written as a switched DAE The coupling via the boundaries of the lines 1, 2 and 3 is as follows: Thus, the overall coupled system has the form (2.1) with where k = 1, 2 and the coupling matrix P = P a 0 0 P b is given by

Linear hyperbolic PDEs
In this section, we recall the solution properties of a system of linear PDEs on a bounded spatial domain where t 0 is the initial time for the system, u : [t 0 , ∞) × [a, b] → R n is the ndimensional unknown vector and A ∈ R n×n with the prescribed initial condition

Assumption 1
The system (3.1) is assumed to be hyperbolic, i.e. A is assumed to be nonsingular and diagonalizable by a real coordinate transformation.
Under Assumption 1, there exists a nonsingular matrix R ∈ R n×n such that where λ i , i = 1, . . . , n, are the (real) eigenvalues of A. The matrix R is composed of the eigenvectors r i corresponding to the eigenvalues˘i of A. Without restricting generality (and under the nonsingularity assumption of A), we can assume that λ 1 ≤ λ 2 ≤ . . . ≤ λ n−r < 0 < λ n−r +1 ≤ . . . ≤ λ n−1 ≤ λ n , where r ∈ {0, 1, . . . , n} is the number of positive eigenvalues of A. Finally, we let := diag(λ 1 , λ 2 , . . . , λ n ) and R = [R − , R + ], where the columns of the matrices R − and R + consist of the eigenvectors that correspond to negative and positive eigenvalues of A, respectively. By applying the coordinate transformation v = R −1 u, we see that (3.1) is equivalent to a decoupled system of scalar PDEs, i = 1, 2, . . . , n, where v t 0 i = r i u t 0 , where r i is the i-th row vector of R −1 . On a bounded spatial domain, say x ∈ [a, b] , it is necessary to prescribe some boundary conditions at the boundaries x = a and x = b. The system (3.1) needs as many boundary conditions at x = a as the number of positive eigenvalues and as many boundary conditions at x = b as the number of negative eigenvalues. The boundary conditions are defined as: where P a ∈ R r ×n , P b ∈ R (n−r )×n , b a : R → R r and b b : R → R n−r . For the wellposedness of the boundary condition (2.1c) for the hyperbolic PDE (3.1) satisfying Assumption 1, we will impose the following assumption on the coupling matrix P.

Assumption 2
The coupling matrix in (2.1c) has the block structure P = P a 0 where P a ∈ R r ×n and P b ∈ R (n−r )×n satisfy ker P a ⊕ im R + = R n and ker For the boundary conditions (3.4) in the characteristic variables, we define M = From Assumption 2, it is easily seen that M 2 and N 1 are invertible, and we therefore can transform (3.4) as follows: v a

Remark 3
In the case that the wave speed λ i = 0, Eq. (3.3) is simply ∂ t v i = 0, meaning that the change in the solution with respect to time is zero, with the initial condition given as in (3.3b). Hence, the solution is given by the initial

Explicit solution formula in terms of characteristic variables
In order to give a solution formula in a compact form, we define two shift operators as follows.

Definition 4 (Shift operator for functions) Denote with F(A → B)
the set of all functions from some set A to some set B. Let T ⊆ R and X ⊆ R, then the time shift operator S λ,x 0 time with speed λ ∈ R and initial position x 0 ∈ X is Consider the system (3.3a) with the initial condition (3.3b) and the boundary conditions (3.5). To ease the notation, let K − := {1, . . . , n−r } and K + := {n−r +1, . . . , n}. The solution to each scalar PDE (3.3a) can be found individually in terms of the initial For left-going waves, the solution is of the form: In a similar fashion, the solution to right-going waves is of the form: where j ∈ K + and x ∈ [a, b]. The solutions v − (t, x) and v + (t, x) together form the solution v(t, x) to the system (3.3a) with the initial condition (3.3b) and the boundary conditions (3.5). Hence, v(t, x) can be written as where e i ∈ R (n−r )×(n−r ) and e j ∈ R r ×r are the i-th and j-th directional unit vectors. The solutions at x = a and x = b then are of the form: (3.9)

Solution framework for the linear hyperbolic system
In this section, the solution u(t, x) to the system (3.1) with the initial and boundary conditions (3.2) and (3.4) will be formulated by using the results from Sect. 3.1 and passing to the representation in the original variables, i.e. u = Rv.

Lemma 5
Consider the PDE (3.1) satisfying Assumption 1 and 2 with some given initial trajectory u t 0 as in (3.2) and boundary conditions b a , b b as in (3.4). Let

11b)
and where it is assumed that Then, every classical solution is given by The proof is a direct consequence of the method of characteristics. Details can be found, for example, in [3].

Remark 6
In addition to the 2D shift operator defined as in Definition 4, we analogously define the 1D time shift operator S τ time with τ ∈ R for functions f : T ⊆ R → R as follows: where F is as in Definition 4 and with the convention that f (s where u a (t) and u b (t) given as in (3.10). Below, we express u ab in a compressed form in terms of the 1D time shift . . , n, are given as in (3.11) and the extensions of initial conditions as boundary conditions for the negative times are adapted as in the proof of Lemma 5. Then, the equality (3.13) follows from Eq. (3.10).

Switched differential-algebraic equations
In this section, we consider switched differential-algebraic equations (swDAEs) of the form: is a piecewise constant switching signal with a locally finite set of jump points and is right continuous, Note that the matrix E ξ is not assumed to be non-singular. For the existence and uniqueness of solutions to (4.1), the following definition of regularity of matrix pairs will be employed. Due to the nonsingularity of E ξ , solutions of (4.1) need to satisfy certain algebraic constraints which in general differ from each mode. Hence, at a switching time t s ∈ R the value w(t − s ) may be inconsistent with the mode after the switch, resulting in a jump in w at t s . In fact, derivatives of jumps, i.e. Dirac impulses, may also occur in response to switches [19]. Therefore, the solution space must be enlarged to distributions so that also jumps and Dirac impulses are allowed in the solutions. To this end, piecewisesmooth distributions are considered as the solution space. Below, we first recall the definition of scalar piecewise-smooth functions and distributions as in [17].
where the functions α i : T → R, i ∈ Z, are (globally) smooth and the set {t i ∈ T | i ∈ Z} is an ordered discrete set, i.e. t i < t i+1 for all i ∈ Z and the intersection with any compact subset of T only contains finitely many points. The set of all piecewise-smooth functions is denoted by Note that a distribution D τ has point support {τ } if and only if there exist d τ ∈ N, τ is the k-th derivative of the Dirac impulse δ τ at τ ∈ T . It is easily seen that the space of piecewise-smooth distributions is closed with respect to differentiation and contains the space of piecewise-smooth functions as a subspace. In fact, a distribution is piecewise-smooth if, and only if, locally it is equal to a finite derivative of a piecewisesmooth function; in other words, for all D ∈ D pwC ∞ and all compact subsets K ⊆ T there exists k ∈ N and α ∈ C ∞ pw (T ; R) such that Having defined a suitable solution space D pwC ∞ that allows distributions as solutions to the swDAE (4.1), the state variable w, the input q and the inhomogeneity f in the swDAE (4.1) are considered in this solution space D pwC ∞ and they are vectors of distributions in D pwC ∞ , i.e. w, f ∈ D pwC ∞ m and q ∈ D pwC ∞ ν .
Theorem 10 (Existence and uniqueness of solutions of swDAEs, [17]) Consider the swDAE as given in (4.1) with regular matrix pairs (E ξ , H ξ ) with ξ ∈ {1, 2, . . . , N }, and assume that the switching signal σ only has locally finitely many switches. Then, for every initial trajectory w t 0 ∈ D pwC ∞ m with the initial time t 0 ∈ R, any input q ∈ D pwC ∞ ν and any inhomogeneity f σ ∈ D pwC ∞ m , there exists a unique globally defined solution w ∈ D pwC ∞ m of the initial trajectory problem Furthermore, the solution on [t 0 , ∞) is uniquely determined by w(t − 0 ).

Distributional solution of the PDE
In Sect. 3, we have reviewed the classical solution to linear hyperbolic PDEs. But considering a coupling with switched DAEs, the boundary data for the PDE are given by piecewise-smooth distributions. Thus, we need to extend the solutions in the distributional sense, including Dirac impulses and its derivatives. Unfortunately, we cannot simply consider distributions on R 2 , since we still need to evaluate the traces at initial time and the boundaries. Therefore, we construct an appropriate solution space by piecewise-smooth distributions in time and space.

Distribution theory in time and space
Definition 11 (Piecewise-smooth functions in time and space) Denote by T ⊆ R (time) and X ⊆ R (space) open intervals. We say a family of subsets (P i ) i∈I of T × X for some index set I is a polyhedral partition of T × X if and only if P i are polyhedral sets (i.e. the intersection of finitely many open or closed half-spaces in T × X ), are pairwise disjoint and i∈I P i = T × X . A function β : T × X → R is called (polyhedral) piecewise-smooth if and only if there exists a locally finite polyhedral partition i∈I P i of T × X and a family of smooth functions β i : where χ P i is the characteristic function of the set P i ⊆ T × X .
In the definition of a polyhedral partition (P i ) i∈I , it is not excluded that some P i have empty interior (i.e. have measure zero), and this has the advantage that then the corresponding space of piecewise-smooth functions is closed under addition and restriction to polyhedral sets. For a piecewise-smooth function β : T × X → R, it is easily seen that for any t ∈ T and x ∈ X , the functions β(t, ·) and β(·, x) are scalar piecewise-smooth functions as in Definition 9 (where we treat two functions as equal when they are equal almost everywhere).
Definition 12 (Dirac segment, cf. [23]) Let L ⊆ T × X be a line segment, i.e. there exists t 0 , t 1 ∈ T , x 0 , x 1 ∈ X such that Then, the Dirac segment on L is where L ϕ is the usual line integral given by where t = t 1 − t 0 and x = x 1 − x 0 . For unbounded line segments (i.e. ξ ranges over an unbounded interval in (5.2)), the integral boundaries in the definition of L ϕ are replaced by ±∞ appropriately.
Note that if t = 0, then and if x = 0 then Lemma 13 Assume T = R, X = R and consider the unbounded line L = {(t 0 + λ t, x 0 + λ x) | λ ∈ R} for some t 0 ∈ T , x 0 ∈ X and t > 0, x > 0. For the step function along L given by Recall the general definition of the partial derivative of a distribution D on T × X : Hence, we have From this, the claims follow.
The space of piecewise-smooth distributions on T × X is denoted by D pwC ∞ (T × X ) Lemma 16 Let D ∈ D pwC ∞ (T × X ) given by (5.3) and (5.1), then The restriction of D to any polyhedral set P ⊆ T × X given by D P := ( i∈I χ P i ∩P β i ) D + j∈J k, α k, j (∂ t ) k (∂ x ) δ L j ∩P is well defined and is again a piecewise-smooth distribution. Proof 1. It suffices to show that the (partial) derivative of (the induced distribution by) a piecewise-smooth function (in the sense of Definition 11) is a piecewisesmooth distribution. Since the sum in (5.1) is locally finite, it furthermore suffices to consider only a single summand and since the multiplication with a smooth function is well defined for general distributions, it suffices to show that the partial derivatives of the indicator function χ P for any polyhedral set P ⊆ T × X are a piecewise-smooth distribution, which was already established in Corollary 14. 2. First note that the intersection of two polyhedral sets is again a polyhedral set; hence, i∈I χ P i ∩P β i is a piecewise-smooth function. Furthermore, the intersection of a line-segment with a polyhedral set is again a line-segment (or empty); hence, D P is again a piecewise-smooth distribution (taking into account the local finiteness property, which implies that evaluated at any test-function reduces to finite sums for which no convergence issues occur due to the restriction).
For any (t, x) ∈ T × X , we want to define in the following the partial evaluations D(t + , ·), D(t − , ·), D(·, x − ), D(·, x + ) such that they are piecewise-smooth distributions on X or T , respectively, and such that (partial) differentiation commutes with the partial evaluation, i.e.
here (·) denotes the (scalar) differentiation in D pwC ∞ (T ) or D pwC ∞ (X ), respectively. Clearly, for piecewise-smooth functions, such an evaluation is trivially defined. Due to commutativity requirement concerning differentiation and evaluation, it is also clear that it suffices to define the evaluation of Dirac segments. Due to Lemma 13, there is, however, only one possible choice:

Distributional solutions for linear hyperbolic PDE
Before addressing linear systems, we consider the scalar advection equation where λ ∈ R is the wave speed and the initial condition is prescribed as where v t 0 ∈ D pwC ∞ ((a, b)) and the boundary condition given as We now expand the definition of the shift operator for continuous functions in Definition 4 to distributions.

Definition 18 (Shift operator for Dirac impulses) Let T , X ⊆ R be open intervals.
The distributional time shift operator of a Dirac impulse at t * ∈ T δ t * ∈ D pwC ∞ (T ) with speed λ and initial position x 0 is given by the distributional space shift operator of a Dirac impulse δ x * ∈ D pwC ∞ (X ) at x * ∈ X with speed λ and initial time t 0 is given by Note that in the definition of the shift operator for Dirac impulses, the factors 1/ 1 + 1/λ 2 and 1/ √ 1 + λ 2 are necessary to obtain the following equalities: S λ,x 0 time δ t * (·, x ± ) = δ t * +(x−x 0 )/λ and S λ,t 0 space δ x * (t ± , ·) = δ x * +λ(t−t 0 ) .
An illustration of the time-and space shift of Dirac impulses can also be found in Fig. 5.

Definition 19 (Shift operator for piecewise-smooth distributions) Let D T ∈ D pwC ∞ (T )
and D X ∈ D pwC ∞ (X ) be given by , T * ⊂ T and X * ⊂ X are locally finite sets, and for each t * ∈ T * and x * ∈ X * we have n t * ∈ N and n x * ∈ N, c t * i ∈ R, i = 0, 1, . . . , n t * and c x * j ∈ R, j = 0, 1, . . . , n x * such that Then, the distributional time shift of D T with speed λ and initial position x 0 is given by and the distributional space shift of D X with speed λ and initial time t 0 is given by Assume λ > 0 and let v t 0 and v a be given as in (5.6) and (5.7), respectively. Below, we will formulate the solution to the Eq. (5.5) in terms of the distributional space/time shift operators S λ,t 0 space , S λ,a time . As seen in Sect. 3, since the solution is constant along the characteristics, exploiting the distributional space/time shift operator given in Definition 19, it can be written as Then, the solution to the differential Eq. (5.5) at the right boundary with λ > 0 takes the form: which can be put in the form a, b). Now assume λ < 0 and let v t 0 and v b be given as in (5.6) and (5.7), respectively, and let v t 0 ∈ D pwC ∞ (X ) and v b ∈ D pwC ∞ (T ). The solution formulae to Eq. (5.5) with λ < 0 are now of the form: Similarly, the solution to the differential Eq. (5.5) with λ < 0 at the left boundary can be written as: which is written as 0 outside (a, b).
As a system of PDEs with boundary conditions, we consider and right boundary conditions with P a ∈ R r ×n , and P b ∈ R (n−r )×n . As in Sect. 3 and with Assumption 1, the system is decomposed into its distributional where the boundary conditions for the right-and left-going waves can be expressed as: a, b)) n−r stands for the left-going waves, whilst v + ∈ D pwC ∞ ((t 0 , ∞) × (a, b)) r for the right-going waves.
The distributional solution v to the decomposed system of (5.11a), as was done similarly to (3.3a), with the initial condition (5.12) and boundary conditions (5.13a)-(5.13b) can be written in terms of the solutions of the left-and right-going waves (5.14) The distributional solution u to the IBVP system (5.11a)-(5.11b)-(5.11c) is now formulated via inversion of the distributional characteristic variables u = Rv. Let 5 p := R diag(e p ) R −1 with diag(e p ) ∈ R n is the p-th directional unit vector. The solution is or, in the compact form with the convention that u t 0 = 0 outside (a, b).
At the left and right end of the spatial domain, the distributional solution u is as follows: where u a := u(·, a + ), u b := u(·, b − ) and with the convention that u t 0 = 0 outside (a, b).

Remark 20
Similar to the 1D time shift defined in Remark 6, we define the 1D distributional time shift S τ time for D ∈ D pwC ∞ (T ). Let D ∈ D pwC ∞ (T ) be given by pw (T ), T * ⊂ T is a locally finite set, and for each t * ∈ T * we have n t * ∈ N, c t * i ∈ R, i = 0, 1, . . . , n t * such that Then, the 1D distributional time shift of D is given by

Remark 21
Equation (5.17) can now be written in compressed form in terms of u ab ∈ D pwC ∞ (T ) 2n and the 1D distributional time shift S τ time as the matrices F a , F b , D ab k , D ba k are given as in (3.11). The extension of the initial conditions as the boundary conditions for the negative times is described in (5.18). Hence, the equality (5.19) follows from the equations in (5.17).
So far we have constructed a piecewise-smooth distributional solution to the hyperbolic PDE (5.11a). This solution is unique due to the following theorem.

Theorem 22 (Uniqueness of the distributional solution)
The solutions to (5.11a) are unique in the space of piecewise-smooth distributions.
This expression is only zero for all ϕ, if λ i = x t and (t 0 , x 0 ) as well as (t 1 , x 1 ) are outside of the support of the ϕ. Thus, the line has to have the slope according to the characteristic speed, and the line has to fully cross the considered domain. But at the points where the line hits the initial time or the boundaries of the domain, it has to satisfy the imposed conditions. Therefore, the strength of the impulse is equal to zero, i.e. the factor c for δ L is zero. Due to linearity, the above computation can be extended directly to any combination of spatial and temporal derivatives of δ L , which concludes the proof.

Existence and uniqueness of the coupled system
In this section, we consider the switched DAE (2.1d) with output y D given by (2.1e) together with the boundary behaviour u ab := u(·, a) , u(·, b) of the PDE (2.1a). Based on the results from the previous section, we can now relate the solution w and u ab of the coupled system one-to-one with the solution of a switched delay DAE as follows.
where τ k , u ab and the matrices F and D k are given as in Remark 21 and its components as in (3.11).
Proof (⇐) Assume that z = (w , u ab ) solves the swDDAE (6.1). From the swD-DAE (6.1), we obtain for which w is the solution with the input q = C P u ab as shown in Theorem 10. And also, from the swDDAE (6.1), we obtain with the convention that u t 0 = 0 outside (a, b) and 5 p := R diag(e p ) R −1 with diag(e p ) ∈ R n is the p-th directional unit vector.
(⇒) Assume that w is the solution to the swDDAE (2.1d) where y P = C P u ab and that u is the solution to the PDE (2.1a) given as (5.15). Then u ab is of the form Hence, z = (w , u ab ) solves the swDDAE (6.1). Since the solution of the coupled system on the whole domain can be recovered via (5.15), we have therefore shown that the solution properties of the coupled system can be equivalently characterized by the swDDAE (6.1). The following result establishes conditions for existence and uniqueness of solutions for general swDDAE.
Theorem 24 (Existence and uniqueness of solutions for swDDAEs) Consider the following switched delay differential-algebraic equations having d ∈ N delays such that 0 < τ 1 < τ 2 < · · · < τ d Proof The result is a simple consequence from the "method of steps", and the details for DDAEs with a single delay d = 1 can be found in [20]. In the following proof, we adapt the prime notation ( ) to indicate the derivative of distributions in addition to the dot notation (˙) since the derivative of a distribution restricted to an interval and restricting a derivative of a distribution to an interval do not have the same meaning. In other words, the operations restriction to an interval and differentiation of a distribution do not commute. Let τ be the smallest delay within the set {τ 1 , . . . , τ d }. The solution to the swDDAE system (6.3) is shown to be expressed as where t k := t 0 + kτ , k ∈ N, t 0 := t 0 , and z k ∈ D pwC ∞ m is the unique solution to the non-delay swDAE, (Theorem 10), For each ϕ ∈ C ∞ pw , the formally infinite sum defining z(ϕ) is actually a finite sum, because ϕ has a compact support. Moreover, z ∈ D pwC ∞ m since it is a linear combination of piecewise-smooth distributions. For any k ≥ 1, where we exploit the relations from cf. [17], as (D [s,t)

Remark 25
The existence and uniqueness result of Theorem 24 can easily be extended to the case that the delay coefficient matrices D j , j = 1, 2, . . . , d in (6.3) are switch dependent, i.e. (6.3) becomes We can conclude now our main result about existence and uniqueness of solutions of the coupled system.
Corollary 26 Consider the coupled system (2.1) with a hyperbolic PDE (Assumption 1) and suitable boundary condition (Assumption 2). Furthermore, assume that for all ξ ∈ {1, . . . , N } the matrix pairs (E ξ , H ξ + B ξ C P FC D ξ ) with F as in Remark 21 are regular. Then, for any initial values Proof This is a consequence of Theorems 23 and 24 and the fact that det( Remark 27 For general delay differential equations, it is common to distinguish the three types retarded, neutral and advanced, see, for example, [9,22]. For our existence and uniqueness result, we do not need any assumption on these different types and in general all three types are allowed. However, when numeric considerations are important, then it is necessary to further analyse the properties of the coupled systems with respect to the index and type of the resulting delay DAE; in particular, the question whether discontinuities in the initial data get smoothed out over time or not is highly relevant. This question was already investigated for nonswitched DDAEs [21] and it should be possible to extend these results to the switched case; this is ongoing research.

Remark 28
As illustrated already in Example 2, even if the PDE is well-defined and the swDAE is regular, coupling these systems may yield an ill-posed system. Indeed, using the notation of Corollary 26, we see that F = 1 0 (note that F b is a 1 × 0 matrix) and hence which is not a regular matrix pair anymore. Furthermore, the equivalent delay DAE according to Theorem 23 does not contain a delay term; hence, the closed-loop systems results in an DAE which is not regular.

Numerical results of the power grid example
In this section, we explain how we solve the coupled PDE system (2.5) with the switched DAE (2.10) numerically and illustrate the results obtained. We remark that the matrix pairs (E σ , H σ ) in the power grid example are regular. Denote by E := {1, 2, 3, 4} the set of edges in the coupled network, over which we discretize the PDE. For discretizing the spatial domain [a, b], we consider x k = a + k x, k = 0, 1, . . . , N , x = (b − a)/N , where N is the number of cells in the mesh. We then insert two ghost cells G −1 and G N +1 at both ends of the computational domain which are treated as boundaries. The discretization of the time variable is done in accordance with the CFL number, which we choose to be 1, t n+1 = t n + t. We denote by u n j = u(x j , t n ) the approximated solution at time t n and position x j . For the power grid example, denote by to have a suitable representation for the unknowns for every edge k ∈ E. Before we start solving the coupled system at each time, we first decompose the discretized PDE over each edge k ∈ E into its characteristic variables (v k ) n j as left-going, (v − k ) n j , and right-going characteristic waves, To solve the decomposed PDE and the swDAE numerically, we use the upwind scheme and implicit Euler method, respectively. At each time iteration, we solve the decomposed PDE numerically, then we update (u k ) n+1 where A k is the coefficient matrix given as in (2.9). The equations in (2.7) together result in coupling conditions in terms of the characteristic variables as follows: which are four out of eight of boundary conditions for the decoupled PDE, and hence, they build up the inputs to the swDAE as The remaining four inputs to the swDAE are given in terms of characteristic vari- Then, at each time step, we solve the swDAE and obtain the boundary conditions Fig. 4, the boundary condition is assigned as (u 2 1 ) where v G is prescribed constant voltage source, and hence, the boundary condition for the characteristics (v + 1 ) n −1 . At P 2 over E 1 , the input to the swDAE is n N +1 and the boundary condition for the characteristic variable is (u 2 1 ) n N +1 . The swDAE assigns this boundary condition where η = 2 or η = 3 depending on to which edge the switch connects E 1 and i n 12 , i n 13 , v n 12 , v n 13 are discretized state variables for the swDAE at time t n . If η = 2, v n+1 12 is computed as above, then the boundary condition (v − 1 ) n N +1 = (v + 1 ) n N − v n+1 12 is assigned; hence, i n+1 12 = (v + 1 ) n N + (v − 1 ) n N +1 . Furthermore, i n+1 13 = 0 and n 0 are assigned. Then, we update the solution for (u k ) n+1 j , for all k ∈ E by using the eigenvector matrix R k and (v − k ) n+1 , = 0, 1, . . . , N and (v + k ) n+1 ψ , ψ = 1, 2, . . . , N + 1. Until we reach the prescribed final time, we update the time t n+1 = t n + t and repeat the above steps. If, instead of solving for (v k ) n j , one attempts to solve for (u k ) n j and assigns inputs/outputs without considering them in characteristic variables, oscillations might occur at boundaries at each time step. The method described in this section covers the Dirac impulses and ensures that such oscillations do not take place. Therefore, the numerical steps defined above should be carried out in characteristic variables, and then, the original unknown variables (u k ) n j should be updated accordingly.

Discontinuous initial condition
We consider the computational domain . For t ∈ [0.7, 1.5), it connects the edges 1 and 2 again. In Fig. 6, the plots for E 1 and E 2 are shown at t = 0.5. After the first switch at t = 0.4, a Dirac impulse occurs on E 2 . In Fig. 7, the plots for all edges at t = 0.8 are shown. After the second switch at t = 0.7, there happens another Dirac impulse on E 3 .

Conclusion
We have presented an existence and uniqueness result for solutions of linear hyperbolic PDE coupled to a linear switched DAE. Since switched DAEs can have distributional solutions (including Dirac impulses and their derivatives of arbitrary high order), it was necessary to first define a suitable distributional solution space for the PDE such that evaluation at the boundary (the trace operator) is well defined. Within this solution space (the space of piecewise-smooth distributions), it is possible to relate the solutions of the coupled systems with solutions of a delay switched DAE and we obtain our main result by extending recent results about the solvability of delay DAEs.
More research is still necessary to extend the results to PDEs with a source term as well as to the mildly nonlinear case. Furthermore, the well-posedness condition we present is still rather technical and there may be simpler sufficient condition, e.g. in terms of the coupling and output matrices of the PDE alone (which then rules out a situation as in Remark 28).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.