Simulation of electromagnetic descriptor models using projectors

Abstract

Electromagnetic descriptor models are models which lead to differential algebraic equations (DAEs). Some of these models mostly arise from electric circuit and power networks. The most frequently used modeling technique in the electric network design is the modified nodal analysis (MNA) which leads to differential algebraic equations in descriptor form. DAEs are known to be very difficult to solve numerically due to the sensitivity of their solutions to perturbations. We use the tractability index to measure this sensitivity since it can be computed numerically. Simulation of DAEs is a very difficult task especially for those with index greater than one. To solve higher-index DAEs, one needs to use multistep methods such as Backward difference formulas (BDFs). In this paper, we present an easier method of solving DAEs numerically using special projectors. This is done by first splitting the DAE system into differential and algebraic parts. We then use the existing numerical integration methods to approximate the solutions of the differential part and the solutions of the algebraic parts are computed explicitly. The desired solution of the DAE system is obtained by taking the linear combination of the solutions of the differential and algebraic parts. Our method is robust and efficient, and can be used on both small and very large systems.

1 Introduction

Consider a linear Resistor-Inductor-Capacitor (RLC) electric network which connects linear capacitors, inductors and resistors, and independent voltage $\mathit{v}(t)\in {\mathbb{R}}^{n_V}$ and current sources $ı(t)\in {\mathbb{R}}^{n_I}$. Circuit designers are interested to know the unknowns which describe the network which are: the node potentials $e(t)\in {\mathbb{R}}^{n_e}$, and the currents ${ı}_L(t)\in {\mathbb{R}}^{n_L}$ and ${ı}_V(t)\in {\mathbb{R}}^{n_V}$ through inductors and voltage sources, respectively. Using the commonly used Modified Nodal Analysis (MNA) [1], we introduce the incidence matrices ${\mathbf{A}}_C\in {\mathbb{R}}^{n_e,n_C}$, ${\mathbf{A}}_L\in {\mathbb{R}}^{n_e,n_L}$ and ${\mathbf{A}}_R\in {\mathbb{R}}^{n_e,n_G}$, which describe the branch-node relationships for capacitors, inductors and resistors. Furthermore, we introduce the incidence matrices ${\mathbf{A}}_V\in {\mathbb{R}}^{n,n_V}$ and ${\mathbf{A}}_I\in {\mathbb{R}}^{n_e,n_I}$, which describe this relationship for voltage and current sources, respectively. This leads to a linear time invariant (LTI) DAE system in descriptor form:

$$\underset{\mathbf{E}}{\underbrace{\left[\begin{array}{ccc}{\mathbf{A}}_C\mathbf{C}{\mathbf{A}}_C^T& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{L}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\end{array}\right]}}{x}^{\prime }=\underset{\mathbf{A}}{\underbrace{\left[\begin{array}{ccc}-{\mathbf{A}}_R\mathbf{G}{\mathbf{A}}_R^T& -{\mathbf{A}}_L& -{\mathbf{A}}_V\\ {\mathbf{A}}_L^T& \mathbf{0}& \mathbf{0}\\ {\mathbf{A}}_V^T& \mathbf{0}& \mathbf{0}\end{array}\right]}}x+\underset{\mathbf{B}}{\underbrace{\left[\begin{array}{cc}-{\mathbf{A}}_I& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\\ \mathbf{0}& -\mathbf{I}\end{array}\right]}}\underset{u}{\underbrace{\left[\begin{array}{c}ı\\ \mathit{v}\end{array}\right]}},$$

(1)

where $\mathbf{C}\in {\mathbb{R}}^{n_C,n_C}$, $\mathbf{L}\in {\mathbb{R}}^{n_L,n_L}$ and $\mathbf{G}\in {\mathbb{R}}^{n_G,n_G}$ are the capacitance, inductance and conductance matrices, respectively which are usually assumed to be symmetric and positive-definite. We need to find the unknowns $x={(e,{ı}_L,{ı}_V)}^T$ of the system (1). If we let $n=n_e+n_L+n_V$ and $m=n_I+n_V$, then $\mathbf{E},\mathbf{A}\in {\mathbb{R}}^{n,n}$, $\mathbf{B}\in {\mathbb{R}}^{n,m}$ and $u\in {\mathbb{R}}^m$. Equation (1) can be solved depending to the elements in the electric networks. For example if we consider a resistor network, (1) simplifies to a system of algebraic equations which implies that $\mathbf{E}=\mathbf{0}$. Thus we need to solve a linear system of the form $-\mathbf{A}x=\mathbf{B}u$. But if the network contains a combination of resistors and other elements, this leads to differential algebraic equations (DAEs), i.e., $det(\mathbf{E})=0$. In this paper, we consider electric networks which leads to DAEs.

DAEs are very difficult to solve numerically due to the sensitivity of their solutions to perturbations. This sensitivity is measured by the index concepts such as differentiation index, perturbation index and tractability index. We use the tractability index because it can be obtained numerically. According to [1], a DAE-index from electric networks cannot be greater than 2. Thus we assume (1) has maximum index of 2. In order to solve DAEs we assume the following: (i) the matrix pencil $\lambda \mathbf{E}-\mathbf{A}$ must be nonsingular for some $\lambda \in \mathbb{C}$. (ii) the input vector u must be smooth enough. Further more consistent initial values $x(0)=x_0$ must be applied.

A lot of work has been done to solve DAE system using multistep methods such BDFs, NDFs and Runge–Kutta methods [2]. Although these methods are accurate they need a lot of computational effort depending on the index of DAE system, for example index-2 systems, BDF is convergent and globally accurate to $\mathcal{O}(h^k)$ but these require tight Newton solutions accurate to $\mathcal{O}(h^{k+1})$ [3]. The implicit Euler method loses order of accuracy with each increase in index, thus cannot be used in practice to solve DAE systems numerically [3].

In this paper, we propose a new way of solving DAEs using special projectors. This method enables one to use any ODE numerical methods to solve DAE systems. This idea was first introduced by März [4], and it involves splitting the DAE system into differential and algebraic parts using projector and matrix chains. Then one can use existing numerical integrations methods to solve the differential part and the some components of the algebraic part are approximated using numerical differential methods for DAEs with index greater than one. The desired solution is obtained by the linear combination of the solutions of the differential and algebraic parts. This approach leads to a decoupled system of dimension $(\mu +1)n$, where μ is the tractability index. Moreover, the spectrum of the decoupled system consists not only of the spectrum of the matrix pencil $(E,A)$ of the original system but also of additional infinite eigenvalues [5]. In practice a system should be stable otherwise it may not work properly or may even be destroyed in practical use. There is a direct criterion for the stability of regular descriptor linear systems which indicates that the stability of the system (1) is totally determined by spectrum of the matrix pencil $\sigma (\mathbf{E},\mathbf{A})$ [6]. The set of all finite eigenvalues of the matrix pencil $(\mathbf{E},\mathbf{A})$ is denoted by ${\sigma }_f(\mathbf{E},\mathbf{A})$ which is called the finite spectrum while the infinite spectrum of the matrix pencil $(\mathbf{E},\mathbf{A})$ is denoted by ${\sigma }_{\mathrm{\infty }}(\mathbf{E},\mathbf{A})$. Thus the spectrum of $(\mathbf{E},\mathbf{A})$ is donated by $\sigma (\mathbf{E},\mathbf{A})={\sigma }_f(\mathbf{E},\mathbf{A})\cup {\sigma }_{\mathrm{\infty }}(\mathbf{E},\mathbf{A})$. System (1) is stable if and only if

$${\sigma }_f(\mathbf{E},\mathbf{A})\subset {\mathbb{C}}^{-}=\{s\in \mathbb{C},Re(s)\le 0,\text{ and any }Re(s)=0\text{ is simple}\}.$$

Thus, we need to make sure the physical properties of the DAE system (1) such as stability to be inherited in its decoupled system. This motivated us to do some modifications in the März decomposition, using special basis vectors, which leads to a modified decoupled system of dimension n. Moreover, this decoupling preserves the spectrum of the matrix pencil $(E,A)$ of the DAE.

This paper is organized as follows: In Section 2, we discuss the decoupling of LTI DAE system using März decomposition and its limitation. In Section 3, we propose the modification of the März decomposition using bases of special projectors. In this section we further discuss the modified decoupling of index-1 and -2 electric systems. We observed that higher index DAEs can decoupled into two ways depending on the spectrum of the matrix pencil $(\mathbf{E},\mathbf{A})$. In Section 4, we check whether the decoupled system preserves the physical properties of the DAE system. In Section 4.1, we compare the numerical accuracy of the decoupled and undecoupled system. In Section 5, we test the proposed method on both simple and industrial problems. The industrial examples show the feasibility of this method on real-life applications. This paper is concluded by some final remarks in Section 6.

2 Decoupling of LTI DAE systems using projectors

Assume (1) is of tractability index μ. Setting ${\mathbf{E}}_0:=\mathbf{E}$, ${\mathbf{A}}_0:=\mathbf{A}$, then the projector and matrix chains of the matrix pair $(\mathbf{E},\mathbf{A})$ can be written as: ${\mathbf{E}}_{j+1}={\mathbf{E}}_j-{\mathbf{A}}_j{\mathbf{Q}}_j$, ${\mathbf{A}}_{j+1}:={\mathbf{A}}_j{\mathbf{P}}_j$, $j\ge 0$, where the image space of the projector ${\mathbf{Q}}_j$ is equal to the nullspace of singular matrix ${\mathbf{E}}_j$, i.e., $Im{\mathbf{Q}}_j=Ker{\mathbf{E}}_j$ and ${\mathbf{P}}_j=\mathbf{I}-{\mathbf{Q}}_j$ is the complementary projector. There exists μ such that ${\mathbf{E}}_{\mu }$ is nonsingular while ${\mathbf{E}}_j$ are singular for all $0\le j<\mu -1$. In order to make sure projector products are also projectors, März [4] suggested an additional constraint ${\mathbf{Q}}_j{\mathbf{Q}}_i=0$, $j>i$ on the projector construction for higher index DAEs. Using these chains we can rewrite (1) as projected system of index-μ:

$${\mathbf{P}}_{\mu -1}\cdots {\mathbf{P}}_0{x}^{\prime }+{\mathbf{Q}}_{0}x+\cdots +{\mathbf{Q}}_{\mu -1}x={\mathbf{E}}_{\mu }^{-1}({\mathbf{A}}_{\mu }x+\mathbf{B}u).$$

(2)

As a consequence, (2) can be decoupled into 1 differential part and μ algebraic parts. Unfortunately, März decomposition leads to decoupled system of larger dimension $(1+\mu )n$ than the dimension n of (1) and also it doesn’t preserve the spectrum of the matrix pencil. This motivated us to modify her decoupled system using special basis vectors introduced in [7] and [5] for index-1 and -2, respectively. This modification leads to a decoupled system which preserves the dimension and the stability of the DAE system. In fact, the dimension of the differential part is equal to the dimension of the finite spectrum of the matrix pencil $(\mathbf{E},\mathbf{A})$. Thus, the stability of the solutions of the decoupled system is guaranteed. We discuss this modification for index-1 and -2 systems in the next section.

3 Modification of the März decomposition

In this section, we modify the März decomposition using special basis vectors. This decomposition leads to a decouple system that preserves the dimension and the spectrum of the matrix pencil of the DAE system.

3.1 Decoupling of index-1 electric networks

Here we assume Equation (1) is of index-1 which implies $\mu =1$ in Equation (2). We then construct basis vectors $({\mathbf{p}}_0,{\mathbf{q}}_0)$ in ${\mathbb{R}}^n$ for the projectors ${\mathbf{P}}_0$ and ${\mathbf{Q}}_0$ with their inversion ${({\mathbf{p}}_0^{\ast },{\mathbf{q}}_0^{\ast })}^T$, where ${\mathbf{p}}_0^T\in {\mathbb{R}}^{n,n_0}$, ${\mathbf{q}}_0^T\in {\mathbb{R}}^{n,k_0}$ and $n=n_0+k_0$. Thus an index-1 electric network can be decoupled as:

(3a)

(3b)

Then, the solution of (1) can be computed using the formula:

$$x={\mathbf{p}}_0{\xi }_p+{\mathbf{q}}_0{\xi }_q,$$

(4)

where ${\mathbf{A}}_p:={\mathbf{p}}_0^{\ast T}{\mathbf{E}}_1^{-1}{\mathbf{A}}_0{\mathbf{p}}_0\in {\mathbb{R}}^{n_0,n_0}$, ${\mathbf{B}}_p:={\mathbf{p}}_0^{\ast T}{\mathbf{E}}_1^{-1}\mathbf{B}\in {\mathbb{R}}^{n_0,m}$, ${\mathbf{A}}_q:={\mathbf{q}}_0^{\ast T}{\mathbf{E}}_1^{-1}{\mathbf{A}}_0{\mathbf{p}}_0\in {\mathbb{R}}^{k_0,n_0}$ and ${\mathbf{B}}_q:={\mathbf{q}}_0^{\ast T}{\mathbf{E}}_1^{-1}\mathbf{B}\in {\mathbb{R}}^{k_0,m}$. Observe that (3a) and (3b) can be solved in a hierarchical way and the desired solutions of the DAE system can be obtained using (4). We observe that the differential part (3a) can be solved using ODE numerical methods and algebraic part (3b) can solved explicitly using the solutions from the differential part. This decoupled system always preserves the dimension since $n=n_0+k_0$ and the stability of (1) since it can easily be proved that $\sigma ({\mathbf{A}}_p)={\sigma }_f(\mathbf{E},\mathbf{A})$. The proof can be found in [7]. $n_0$ and $k_0$ is the number of differential and algebraic equations, respectively.

3.2 Decoupling of index-2 electric networks

We now assume (1) is of index-2 which implies $\mu =2$ in Equation (2). We first construct basis vectors $({\mathbf{p}}_0,{\mathbf{q}}_0)$ in ${\mathbb{R}}^n$ with their inversion ${({\mathbf{p}}_0\ast ,{\mathbf{q}}_0\ast )}^T$ for the projectors ${\mathbf{P}}_0$ and ${\mathbf{Q}}_0$, where ${\mathbf{p}}_0\in {\mathbb{R}}^{n,n_0}$, ${\mathbf{q}}_0\in {\mathbb{R}}^{n,k_0}$. For this case we have two possibilities depending on the spectrum of the matrix pencil $(\mathbf{E},\mathbf{A})$.

3.2.1 Matrix pencil $(\mathbf{E},\mathbf{A})$ with at least one finite eigenvalue

Here, we assume that the matrix pencil $(\mathbf{E},\mathbf{A})$ of Equation (1) has at least one finite eigenvalue. Using projector chains ${\mathbf{Q}}_1,{\mathbf{P}}_1\in {\mathbb{R}}^{n,n}$ leads to a theorem below [5].

Theorem 1 Let ${\mathbf{P}}_{01}={\mathbf{p}}_0^{\ast T}{\mathbf{P}}_1{\mathbf{p}}_0$, ${\mathbf{Q}}_{01}={\mathbf{p}}_0^{\ast T}{\mathbf{Q}}_1{\mathbf{p}}_0$. Then ${\mathbf{P}}_{01},{\mathbf{Q}}_{01}\in {\mathbb{R}}^{n_0,n_0}$ are projectors in ${\mathbb{R}}^{n_0}$ provided the constraint condition ${\mathbf{Q}}_1{\mathbf{Q}}_0=0$ holds.

Next, we construct another basis vectors $({\mathbf{p}}_{01},{\mathbf{q}}_{01})$ in ${\mathbb{R}}^{n_0}$ made of $n_{01}$ independent columns of projector ${\mathbf{P}}_{01}$ and $k_1$ independent columns of its complementary projector ${\mathbf{Q}}_{01}$ such that $n_0=n_{01}+k_1$. Let the inverse of this basis be denoted by ${({\mathbf{p}}_{01}^{\ast },{\mathbf{q}}_{01}^{\ast })}^T$. Then Equation (1) can be decoupled as:

(5a)

(5b)

(5c)

and the desired solution of (1) is given by:

$$x={\mathbf{p}}_0{\mathbf{p}}_{01}{\xi }_p+{\mathbf{p}}_0{\mathbf{q}}_{01}{\xi }_{q,1}+{\mathbf{q}}_0{\xi }_{q,0},$$

(6)

where ${\mathbf{A}}_p:={\mathbf{p}}_{01}^{\ast T}{\mathbf{p}}_0^{\ast T}{\mathbf{E}}_2^{-1}{\mathbf{A}}_2{\mathbf{p}}_0{\mathbf{p}}_{01}\in {\mathbb{R}}^{n_{01},n_{01}}$, ${\mathbf{B}}_p:={\mathbf{p}}_{01}^{\ast T}{\mathbf{p}}_0^{\ast T}{\mathbf{E}}_2^{-1}B\in {\mathbb{R}}^{n_{01},m}$, ${\mathbf{A}}_{q,1}:={\mathbf{q}}_{01}^{\ast T}{\mathbf{p}}_0^{\ast T}{\mathbf{E}}_2^{-1}{\mathbf{A}}_2{\mathbf{p}}_0{\mathbf{p}}_{01}\in {\mathbb{R}}^{k_1,n_{01}}$, ${\mathbf{B}}_{q,1}:={\mathbf{q}}_{01}^{\ast T}{\mathbf{p}}_0^{\ast T}{\mathbf{E}}_2^{-1}B\in {\mathbb{R}}^{k_1,m}$, ${\mathbf{A}}_{q,0}:={\mathbf{q}}_0^{\ast T}{\mathbf{P}}_1{\mathbf{E}}_2^{-1}{\mathbf{A}}_2{\mathbf{p}}_0{\mathbf{p}}_{01}\in {\mathbb{R}}^{k_0,n_{01}}$, ${\mathbf{B}}_{q,0}:={\mathbf{q}}_0^{\ast T}{\mathbf{P}}_1{\mathbf{E}}_2^{-1}B\in {\mathbb{R}}^{k_0,m}$ and ${\mathbf{A}}_{q,01}:={\mathbf{q}}_0^{\ast T}{\mathbf{Q}}_1{\mathbf{p}}_0{\mathbf{q}}_{01}\in {\mathbb{R}}^{k_0,k_1}$. Equations (5a), (5b) and (5c) are of dimension $n_{01}$, $k_1$ and $k_0$, respectively, where $n=n_{01}+k_1+k_0$. System (5a)-(5c) also preserves the dimension and stability of system (1). Since it also be proved that $\sigma ({\mathbf{A}}_p)={\sigma }_f(\mathbf{E},\mathbf{A})$. For the purpose of stable numerical implementation decoupled system (5a)-(5c) and (6) can be rewritten as:

(7a)

(7b)

(7c)

where ${\xi }_q={({\xi }_{q,1},{\xi }_{q,0})}^T\in {\mathbb{R}}^{n_q}$, ${\mathbf{A}}_q={({\mathbf{A}}_{q,1},{\mathbf{A}}_{q,0})}^T\in {\mathbb{R}}^{n_q,n_2}$, ${\mathbf{B}}_q={({\mathbf{B}}_{q,1},{\mathbf{B}}_{q,0})}^T\in {\mathbb{R}}^{n_q,m}$, $u^{(i)}=\frac{{\mathrm{d}}^{i}u}{\mathrm{d}t}\in {\mathbb{R}}^m$ and

(8)

is a strictly lower triangular nilpotent matrix of index 2. Also the decoupled system (7a)-(7c) can be solved in a hierarchical way and the desired solutions of the DAE system can be obtained using (7c).

3.2.2 Matrix pencil $(\mathbf{E},\mathbf{A})$ with no finite eigenvalues

We assume that the matrix pencil $(\mathbf{E},\mathbf{A})$ of Equation (1) has no finite eigenvalues, i.e. $det(\lambda E-A)=c\in \mathbb{R}\setminus \{0\}$, $\mathrm{\forall }\lambda \in \mathbb{C}$. Thus Equation (1) can be decoupled as:

(9a)

(9b)

and the desired solution of (1) is given by:

$$x={\mathbf{p}}_0{\xi }_{q,1}+{\mathbf{q}}_0{\xi }_{q,0},$$

(10)

where ${\mathbf{B}}_{q,1}={\mathbf{p}}_0^{\ast T}{\mathbf{E}}_2^{-1}B\in {\mathbb{R}}^{n_0,m}$, ${\mathbf{B}}_{q,0}={\mathbf{q}}_0^{\ast T}{\mathbf{P}}_1{\mathbf{E}}_2^{-1}B\in {\mathbb{R}}^{k_0,m}$ and ${\mathbf{A}}_{q,01}={\mathbf{q}}_0^{\ast T}{\mathbf{Q}}_1{\mathbf{p}}_0\in {\mathbb{R}}^{k_0,n_0}$. Here $n=n_0+k_0$. Still system (9a) and (9b) can simplified to:

(11a)

(11b)

where ${\xi }_q={({\xi }_{q,1},{\xi }_{q,0})}^T\in {\mathbb{R}}^n$, ${\mathbf{B}}_q={({\mathbf{B}}_{q,1},{\mathbf{B}}_{q,0})}^T\in {\mathbb{R}}^{n,m}$, $u^{(i)}=\frac{{\mathrm{d}}^{i}u}{\mathrm{d}t}\in {\mathbb{R}}^m$ and is also a strictly lower triangular nilpotent matrix of index 2 which is defined as (8). We can observe that we do not need initial conditions to solve system (11a) and (11b), in fact the solutions are computed explicitly. We note that the input vector u must be at least $\mu -1$ times differentiable.

We have discussed that index-1 and -2 electric networks can be decoupled using special basis vectors. These special basis vectors can be computed numerically in an efficient way using the LUQ routine discussed in [8]. Hence this procedure can be applied even on very large electric networks as illustrated in Section 5.2. We call this method the Split-DAE method since it involves splitting the DAE system into differential and algebraic variables.

4 Analysis of decoupled electric networks

In this section, we analyse the solutions of the decoupled systems discussed in the previous section. For simplicity we need to first rewrite the decoupled systems back into the descriptor form:

(12a)

(12b)

where $\tilde{\mathbf{E}},\tilde{\mathbf{A}},\mathbf{V}\in {\mathbb{R}}^{n,n}$, $\tilde{\mathbf{B}}\in {\mathbb{R}}^{n,m}$ and $\xi \in {\mathbb{R}}^n$ is the projected state space. We note that the solutions of (1) and (12a) and (12b) coincide since the two systems are equivalent to each other. This implies the spectrum of matrix pair $(\mathbf{E},\mathbf{A})$ and $(\tilde{\mathbf{E}},\tilde{\mathbf{A}})$ must also coincide. System (12a) and (12b) is much easier to solve than its counter part and moreover it reveals the interconnection or structure of the DAE system. Assume the matrix pencil $(\mathbf{E},\mathbf{A})$ of (1) has at least one finite eigenvalue then following theorem below holds:

Theorem 2 If the matrix pair $(\mathbf{E},\mathbf{A})$ has at least one finite eigenvalue, there exist finite $\lambda \in \mathbb{C}$ such that:

$$det(\lambda \mathbf{I}-{\mathbf{A}}_p)=0.$$

This implies that ${\sigma }_f(\mathbf{E},\mathbf{A})=\sigma ({\mathbf{A}}_p)$.

Thus the stability of the decoupled system (12a) and (12b) depends on the stability of the differential part. Hence if Equation (1) is stable then also (12a) and (12b) must be stable.

Using the decoupled systems in Section 3, we derive the matrices $(\tilde{\mathbf{E}},\tilde{\mathbf{A}},\tilde{\mathbf{B}},\mathbf{V})$ below:

• 1. Index 1 electric networks

From decoupled system (3a) and (3b) we obtain:

$$\begin{array}{r}\tilde{\mathbf{E}}:=\left[\begin{array}{cc}{\mathbf{I}}_{n_0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right],\tilde{\mathbf{A}}:=\left[\begin{array}{cc}{\mathbf{A}}_p& \mathbf{0}\\ {\mathbf{A}}_q& -{\mathbf{I}}_{k_0}\end{array}\right],\tilde{\mathbf{B}}:=\left[\begin{array}{c}{\mathbf{B}}_p\\ {\mathbf{B}}_q\end{array}\right],\\ \mathbf{V}:=\left[\begin{array}{cc}{\mathbf{p}}_0& {\mathbf{q}}_0\end{array}\right]\text{and}\xi :=\left[\begin{array}{c}{\xi }_p\\ {\xi }_q\end{array}\right].\end{array}$$

(13)

• 2. Index 2 electric networks

• (i) Matrix pencil $(\mathbf{E},\mathbf{A})$ with at least one finite eigenvalue

From decoupled system (5a)-(5c) we obtain:

(14)

• (ii) Matrix pencil $(\mathbf{E},\mathbf{A})$ with no finite eigenvalues

From decoupled system (9a) and (9b) we obtain:

(15)

We observe that the descriptor form with differential part takes the same form for the case of index-1 and -2 system though for index-1 systems. We use this form only for analysis of the solutions of the DAE system (1) but not for solving. For solving one need to solve the decoupled systems derived in the previous section.

4.1 Numerical accuracy of the decoupled system

Suppose that we want to approximate the solution of the DAE system (12a) and (12b). We can do this by discretizing the time interval $t\in [0,T]$ into N subintervals using $h=\frac{T-0}{N}$ time-steps and set $t_n=ih$, $i=0,1,\dots ,N-1$. We can then approximate the solutions of DAE system (12a) and (12b) at every point. For simplicity, we consider the implicit Euler method since is the starting point for most high order methods and numerically stable. Applying this method on (12a), we obtain:

$$\tilde{\mathbf{E}}\frac{({\xi }_{i+1}-{\xi }_i)}{h}=\tilde{\mathbf{A}}{\xi }_{i+1}+\tilde{\mathbf{B}}{u}_{i+1}.$$

(16)

Rearranging Equation (16) we obtain:

$${\xi }_{i+1}={(\tilde{\mathbf{E}}-h\tilde{\mathbf{A}})}^{-1}[\tilde{\mathbf{E}}{\xi }_i+h\tilde{B}{u}_{i+1}].$$

(17)

Thus the approximate solution of Equations (12a) and (12b) can be obtained using the formula: $x_{i+1}=\mathbf{V}{\xi }_{i+1}$. Then the exact solution of Equations (12a) and (12b) is given by $x(t_i)=\mathbf{V}\xi (t_i)$. The exact solution $\xi (t_i)$ of Equation (12a) satisfies:

$$\tilde{\mathbf{E}}[\frac{\xi (t_{i+1})-\xi (t_i)}{h}+\frac{h}{2}{\xi }^{\mathrm{\prime }\mathrm{\prime }}(\eta )]=\tilde{A}\xi (t_{i+1})+\tilde{\mathbf{B}}u(t_{i+1}).$$

Defining the global error $e_i:=\xi (t_i)-{\xi }_i$ and setting $e_0=0$, since the consistent initial condition is known. We obtain an iterative process below:

$$e_{i+1}={(\tilde{\mathbf{E}}-h\tilde{\mathbf{A}})}^{-1}\tilde{\mathbf{E}}[e_i-\frac{h^2}{2}{\xi }^{\mathrm{\prime }\mathrm{\prime }}(\eta )],i\ge 0.$$

(18)

We use Equation (18) to analyse the accuracy of the implicit Euler method for the first step, i.e. $i=1$ using the following cases.

• 1. Index 1 systems

Using system matrices (13), we obtain:

$${(\tilde{\mathbf{E}}-h\tilde{\mathbf{A}})}^{-1}\tilde{\mathbf{E}}=\left[\begin{array}{cc}{({\mathbf{I}}_{n_0}-h{\mathbf{A}}_p)}^{-1}& \mathbf{0}\\ {\mathbf{A}}_q{({\mathbf{I}}_{n_0}-h{\mathbf{A}}_p)}^{-1}& \mathbf{0}\end{array}\right],{\xi }^{\mathrm{\prime }\mathrm{\prime }}(\eta )=\left[\begin{array}{c}{\xi }_p^{\mathrm{\prime }\mathrm{\prime }}(\eta )\\ {\mathbf{A}}_q{\xi }_p^{\mathrm{\prime }\mathrm{\prime }}(\eta )+{\mathbf{B}}_q{u}^{\mathrm{\prime }\mathrm{\prime }}(\eta )\end{array}\right].$$

(19)

Substituting Equation (19) into Equation (18), leads to:

$$e_1=\frac{h^2}{2}\left[\begin{array}{c}{({\mathbf{I}}_{n_0}-h{\mathbf{A}}_p)}^{-1}{\xi }_p^{\mathrm{\prime }\mathrm{\prime }}(\eta )\\ {({\mathbf{I}}_{n_0}-h{\mathbf{A}}_p)}^{-1}{\xi }_p^{\mathrm{\prime }\mathrm{\prime }}(\eta )\end{array}\right]=\mathcal{O}\left(h^2\right).$$

We observe that this corresponds to the local error for implicit Euler method. This implies that the accuracy of the solutions of the decoupled system is determined by the numerical accuracy of the differential part. As illustrated in Example 1.

• 2. Index 2 systems

As we have already discussed for the case of index-2 system in Section 3.2, we need to consider two cases: We observe that also this case the implicit Euler method loses local error accuracy. For our method we do not need to use any integration method since the solutions of the decoupled system (11a) are computed ellipticity. This is also illustrated in the Example 3.

• (i) Matrix pencil $(\mathbf{E},\mathbf{A})$ with at least one finite eigenvalue

Here we use system matrices (14) and obtain:

(20)

where ${\mathbf{M}}_p={({\mathbf{I}}_{n_{01}}-h{\mathbf{A}}_p)}^{-1}$, and ${\xi }_q(\eta )$ is defined by expression (7b). Substituting Equation (20) into (18), leads to:

(21)

since $h^3{\mathbf{M}}_L{\mathbf{A}}_q{\mathbf{M}}_p=\mathcal{O}(h^2)$ and . We observe that for this case implicit Euler method loses accuracy. We observe that the hidden constraint (7b) destroys the accuracy of the implicit Euler method. We do not face this problem in our method since the algebraic part (7b) is computed explicitly. For our method the accuracy will depend on the method used to solve the differential part (7a). This is illustrated in Example 2.

• (ii) Matrix pencil $(\mathbf{E},\mathbf{A})$ with no finite eigenvalue

Using system matrices (15), we obtain:

(22)

where ${\mathbf{M}}_L$ is defined as in the previous case while ${\xi }_q(\eta )$ is defined by (11a). Substituting Equation (22) into Equation (18), we obtain:

Hence implicit Euler method loses accuracy for each increase in the index if used on the DAE system directly. Thus first decoupling the DAE system (1) into differential and algebraic parts is a robust and effective way of solving DAEs. Since on the differential part one can use any numerical integration method and then the solutions of the algebraic part can be obtained explicitly. Hence the Split-DAE method is a very accurate method and very easy to implement as compared to its counterparts discussed in [2].

5 Numerical experiments

In this section, we test our proposed method in Section 3 using problems from electromagnetic community since is the main focus of this paper but it can also be applied to other applications.

5.1 Simple examples

Here, we explicitly discuss how to decouple DAEs into differential and algebraic parts using projectors. We illustrate this using examples below. Example 1 is an index-1 electric network while Example 2 is an index-2 electric network whose decoupled system has a differential part. In Example 3, we illustrate the decoupling of index-2 system without differential part.

Example 1 Consider a linear RLC electric network in Figure 1. We need to find the unknowns $x={[e_1{e}_2{e}_3{ı}_L{ı}_V]}^T$ in the above electric network. In order to find the unknowns we use the MNA to formulate a mathematical model of the form (1) with system matrices:

$$\mathbf{E}=\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& C& 0& 0\\ 0& 0& 0& L& 0\\ 0& 0& 0& 0& 0\end{array}\right],\mathbf{A}=\left[\begin{array}{ccccc}-G& G& 0& 0& 1\\ G& -G& 0& -1& 0\\ 0& 0& 0& 1& 0\\ 0& 1& -1& 0& 0\\ -1& 0& 0& 0& 0\end{array}\right],\mathbf{B}=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ -1\end{array}\right].$$

(23)

We observe that $det(\mathbf{E})=0$, thus the system is a DAEs of dimension $n=5$. It is solvable since $det(\lambda E-A)=GCL{\lambda }^2+C\lambda +G\ne 0$ for some $\lambda \in \mathbb{C}$. We need to first check the index of the DAE system. We do this by constructing the matrix and projector chains of matrix pair $(\mathbf{E},\mathbf{A})$. Setting ${\mathbf{E}}_0=\mathbf{E}$, ${\mathbf{A}}_0=\mathbf{A}$, we can choose projector ${\mathbf{Q}}_0$ such that $Im{\mathbf{Q}}_0=Ker{\mathbf{E}}_0$, and then compute ${\mathbf{E}}_1={\mathbf{E}}_0-{\mathbf{A}}_0{\mathbf{Q}}_0$. Thus we obtain the first iterate of matrix and projector chain given by:

$$\begin{array}{r}{\mathbf{Q}}_0=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1\end{array}\right],{\mathbf{P}}_0=\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0\end{array}\right],\\ {\mathbf{E}}_1=\left[\begin{array}{ccccc}G& -G& 0& 0& -1\\ -G& G& 0& 0& 0\\ 0& 0& C& 0& 0\\ 0& -1& 0& L& 0\\ 1& 0& 0& 0& 0\end{array}\right].\end{array}$$

(24)

Since ${\mathbf{E}}_1$ is nonsingular, thus Equation (23) is an index-1 system. Thus the decouple system of (23) will take the form (3a) and (3b). We then construct basis matrices and their respective inverses for projectors ${\mathbf{Q}}_0$ and ${\mathbf{P}}_0$ given by:

$${\mathbf{q}}_0=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right],{\mathbf{p}}_0=\left[\begin{array}{cc}0& 0\\ 0& 0\\ 1& 0\\ 0& 1\\ 0& 0\end{array}\right],\text{and}{\mathbf{q}}_0^{\ast T}={\mathbf{q}}_0^T,{\mathbf{p}}_0^{\ast T}={\mathbf{p}}_0^T.$$

(25)

Substituting Equations (23)-(25) into Equations (3a) and (3b) leads to the decoupled system of system (23) given by:

(26a)

(26b)

and the desired solution is given by:

$$x={\left[\begin{array}{ccccc}0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\end{array}\right]}^T{\xi }_p+{\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 1\end{array}\right]}^T{\xi }_q.$$

The decoupled system (26a) and (26b) can be solved in hierarchical way if we apply ${\xi }_p(0)={[e_3(0){ı}_L(0)]}^T$. Observe that $n_0=2$, $k_0=3\Rightarrow n=n_0+k_0=5$ and ${\sigma }_f(\mathbf{E},\mathbf{A})=\sigma ({\mathbf{A}}_p)$ as expected. Thus system (23) is decoupled into 2 and 3 differential and algebraic equations, respectively. We observe that decoupled system (26a) and (26b) can be solved using any ODE numerical integration methods. Thus system (26a) and (26b) is easier to solve than solving (23).

Figure 1

Simple RLC electric network.

Example 2

Consider another linear RLC electric network with system matrices:

$$\begin{array}{r}\mathbf{E}=\left[\begin{array}{cccc}C& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& L& 0\\ 0& 0& 0& 0\end{array}\right],\mathbf{A}=\left[\begin{array}{cccc}-G& G& 0& 1\\ G& -G& -1& 0\\ 0& 1& 0& 0\\ -1& 0& 0& 0\end{array}\right],\\ \mathbf{B}=\left[\begin{array}{c}0\\ 0\\ 0\\ -1\end{array}\right],x=\left[\begin{array}{c}{e}_1\\ e_2\\ {ı}_L\\ {ı}_V\end{array}\right].\end{array}$$

(27)

This is also DAE system of dimension $n=4$ and it is also solvable since $det(\lambda E-A)=GL\lambda +1\ne 0$ for some $\lambda \in \mathbb{C}$. We construct matrix and projector chains by first setting ${\mathbf{E}}_0=\mathbf{E}$, ${\mathbf{A}}_0=\mathbf{A}$. We choose projector ${\mathbf{Q}}_0$ such that $Im{\mathbf{Q}}_0=Ker{\mathbf{E}}_0$ and its complementary projector ${\mathbf{P}}_0$:

$${\mathbf{Q}}_0=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1\end{array}\right],\text{and}{\mathbf{P}}_0=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\end{array}\right].$$

(28)

Then,

$${\mathbf{E}}_1={\mathbf{E}}_0-{\mathbf{A}}_0{\mathbf{Q}}_0=\left[\begin{array}{cccc}C& -G& 0& -1\\ 0& G& 0& 0\\ 0& -1& L& 0\\ 0& 0& 0& 0\end{array}\right],{\mathbf{A}}_1={\mathbf{A}}_0{\mathbf{P}}_0=\left[\begin{array}{cccc}0& G& 0& 1\\ 0& -G& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\end{array}\right].$$

Since ${\mathbf{E}}_1$ is singular, the index-1 condition is violet. Thus we need to continue with the process and choose projector ${\mathbf{Q}}_1$ such that $Im{\mathbf{Q}}_1=Ker{\mathbf{E}}_1$ and its complementary projector ${\mathbf{P}}_1=\mathbf{I}-{\mathbf{Q}}_1$:

$${\mathbf{Q}}_1=\frac{1}{C^2+1}\left[\begin{array}{cccc}1& 0& 0& C\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ C& 0& 0& C^2\end{array}\right],\text{and}{\mathbf{P}}_1=\mathbf{I}-{\mathbf{Q}}_1.$$

(29)

Then,

$${\mathbf{E}}_2={\mathbf{E}}_1-{\mathbf{A}}_1{\mathbf{Q}}_1=\left[\begin{array}{cccc}C+\frac{G}{C^2+1}& -G& 0& -1+\frac{GC}{C^2+1}\\ -\frac{G}{C^2+1}& G& 0& -\frac{GC}{C^2+1}\\ 0& -1& L& 0\\ \frac{1}{C^2+1}& 0& 0& -\frac{C}{C^2+1}\end{array}\right].$$

(30)

We can easily check that ${\mathbf{E}}_2$ is non-singular. Thus Equation (27) is of index-2.

The next step is to decouple the DAE system (27) using the constructed matrix and projector chains. In this step we need to make sure that condition ${\mathbf{Q}}_1{\mathbf{Q}}_0=0$ holds true. We notice that projectors (28) and (29) don’t satisfy this condition. Thus we need to construct other special projectors. Our best choice is the canonical projector suggested by März [4], which is defined as:

$${\tilde{\mathbf{Q}}}_1:=-{\mathbf{Q}}_1{\mathbf{E}}_2^{-1}{\mathbf{A}}_1=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ C& 0& 0& 0\end{array}\right],\text{and}{\tilde{\mathbf{P}}}_1=\mathbf{I}-{\tilde{\mathbf{Q}}}_1.$$

We can easily check that condition ${\tilde{\mathbf{Q}}}_1{\mathbf{Q}}_0=0$ is now satisfied. We then use ${\tilde{\mathbf{Q}}}_1$ instead of ${\mathbf{Q}}_1$ in the decoupling process. If we repeat step (30), we obtain:

$${\mathbf{E}}_2={\mathbf{E}}_1-{\mathbf{A}}_1{\tilde{\mathbf{Q}}}_1=\left[\begin{array}{cccc}C+G& -G& 0& -1\\ -G& G& 0& 0\\ 0& -1& L& 0\\ 1& 0& 0& 0\end{array}\right],{\mathbf{A}}_2={\mathbf{A}}_1{\tilde{\mathbf{P}}}_1=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right].$$

(31)

For convenience we can set ${\mathbf{Q}}_1={\tilde{\mathbf{Q}}}_1$ and ${\mathbf{P}}_1={\tilde{\mathbf{P}}}_1$. Since Equation (27) is an index-2 system and its matrix pencil has at least one finite eigenvalues. Thus its decouple system will be in form (5a)-(5c). The next step is to construct basis matrices for projectors and projector products. The basis matrices for projectors ${\mathbf{P}}_0$ and ${\mathbf{Q}}_0$ in Equation (28), and their respective inverses are given by:

$${\mathbf{p}}_0=\left[\begin{array}{cc}1& 0\\ 0& 0\\ 0& 1\\ 0& 0\end{array}\right],{\mathbf{q}}_0=\left[\begin{array}{cc}0& 0\\ 1& 0\\ 0& 0\\ 0& 1\end{array}\right],\text{and}{\mathbf{p}}_0^{\ast }=\left[\begin{array}{cc}1& 0\\ 0& 0\\ 0& 1\\ 0& 0\end{array}\right],{\mathbf{q}}_0^{\ast }=\left[\begin{array}{cc}0& 0\\ 1& 0\\ 0& 0\\ 0& 1\end{array}\right].$$

(32)

Using Theorem 1, we compute other projectors ${\mathbf{P}}_{01}$ and ${\mathbf{Q}}_{01}$ given by:

$${\mathbf{P}}_{01}={\mathbf{p}}_0^{\ast T}{\mathbf{P}}_1{\mathbf{p}}_0=\left[\begin{array}{cc}0& 0\\ 0& 1\end{array}\right],\text{and}{\mathbf{Q}}_{01}={\mathbf{p}}_0^{\ast T}{\mathbf{Q}}_1{\mathbf{p}}_0=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right].$$

The basis matrices for projectors ${\mathbf{P}}_{01}$ and ${\mathbf{Q}}_{01}$, and their respective inverses are given by:

$${\mathbf{p}}_{01}=\left[\begin{array}{c}0\\ 1\end{array}\right],{\mathbf{q}}_{01}=\left[\begin{array}{c}1\\ 0\end{array}\right],\text{and}{\mathbf{p}}_{01}^{\ast }=\left[\begin{array}{c}0\\ 1\end{array}\right],{\mathbf{q}}_{01}^{\ast }=\left[\begin{array}{c}1\\ 0\end{array}\right].$$

(33)

Substituting Equations (31)-(33) into (7a)-(7c) leads to the decoupled system of (27) given by:

(34a)

(34b)

and the desired solution is given by:

$$x=\left[\begin{array}{c}0\\ 0\\ 1\\ 0\end{array}\right]{\xi }_p+\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right]{\xi }_q.$$

(35)

We observe that $n_{10}=1$, $k_1=1$, $k_0=2\Rightarrow n=n_{10}+k_1+k_0=4$ and ${\sigma }_f(\mathbf{E},\mathbf{A})=\sigma ({\mathbf{A}}_p)=\{-\frac{1}{GL}\}$. Equations (34a) and (34b) can be solved if we apply ${\xi }_p(0)={ı}_L(0)$. Thus the DAE system (27) is decoupled into 1 differential and 3 algebraic equations. The solutions of system (27) and (34a) and (34b), (35) coincides though the later is easier to solve than the former.

Example 3 Consider a simple RL electric network in Figure 2. We need to find the unknowns $x={[e_1,e_2,{ı}_L]}^T$ in the electric network. Using the Modified Nodal Analysis on this network leads to DAE system of the form (1), where $u=ı(t)$ with system matrices given by

$$E=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& L\end{array}\right],A=\left[\begin{array}{ccc}-G& G& 0\\ G& -G& -1\\ 0& 1& 0\end{array}\right],B=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right].$$

(36)

Since $det(\lambda E-A)=G$, thus this system is solvable and its matrix pencil $(\mathbf{E},\mathbf{A})$ has only infinite eigenvalues. Thus its decoupled system must be purely algebraic system. We choose special projectors,

$${\mathbf{Q}}_0=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right],\text{and}{\mathbf{Q}}_1=\left[\begin{array}{ccc}0& 0& L\\ 0& 0& L\\ 0& 0& 1\end{array}\right],$$

(37)

such that ${\mathbf{Q}}_1{\mathbf{Q}}_0=0$ holds. Then we have,

$${\mathbf{E}}_2=\left[\begin{array}{ccc}G& -G& 0\\ -G& G& 1\\ 0& -1& L\end{array}\right].$$

(38)

Since the ${\mathbf{E}}_2$ is non-singular, thus (36) is an index 2 system. The linearly independent columns and their respective inverses of projector ${\mathbf{Q}}_0$ and its complimentary ${\mathbf{P}}_0$ are given by;

$$\begin{array}{r}{\mathbf{p}}_0=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right],{\mathbf{q}}_0=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right],\\ {\mathbf{p}}_0{\ast }^T=\left[\begin{array}{ccc}0& 0& 1\end{array}\right],{\mathbf{q}}_0{\ast }^T=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\end{array}\right].\end{array}$$

(39)

Thus, substituting Equations (37)-(39) into (11a) and (11b), we obtain,

(40)

(41)

We can observe that system (40) can be computed explicitly without using any numerical integration techniques. Thus instead of solving DAE system (36) is easier to solve (40) and their solutions much coincide.

Figure 2

Simple RL network.

5.2 Industrial applications

In this section, we test the Split-DAE method on large scale descriptor electromagnetic models from industries. Examples 4 and 5 are index-1 and -2 systems respectively. These examples are multiple-input multiple-output (MIMO) type, i.e. we are interested in a few solutions of the DAE system.

Example 4 This is a descriptor model of a large scale power system originating from [9]. It is called $\mathrm{b}\mathrm{a}\mathrm{u}\mathrm{r}\mathrm{u}\mathrm{5727}$ power system and can be downloaded from [10]. It is an index-1 DAE system of dimension $n=40\text{,}266$ with 2 inputs and 2 outputs. The sparsity of its matrix pencil is shown in Figure 3. We decoupled the DAE system into $5\text{,}727$ differential equations and $34\text{,}639$ algebraic equations using the procedure we discussed in Section 3.1, i.e. $n_0=5\text{,}727$ and $k_0=34\text{,}639$. The sparsity of the matrix pencil of the decoupled system in the descriptor form is given in Figure 4. We observe that the sparsity in Figure 4 has the same structures as matrix pencil $(\tilde{\mathbf{E}},\tilde{\mathbf{A}})$ in Equation (13). We solved the decoupled system in Matlab using ODE solver ($\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{15}\mathrm{s}$) to solve the differential part. We obtained good solutions as shown in Figure 5 and the simulation took $3\text{,}486$ seconds. We used input function $u(t)={[cos(t)sin(t)]}^T$, $t\in [0,\pi ]$ and time step $h=1.5\times {10}^{-4}$. We tried to solve the DAE system directly using he Matlab DAE solver ($\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{15}\mathrm{s}$) but failed to solve it. It just gave us a warning: ‘This DAE appears to be of index greater than 1’.

Figure 3

Sparsity of matrix pencil $\mathbf{(}\mathbf{E}\mathbf{,}\mathbf{A}\mathbf{)}$ .

Figure 4

Sparsity of matrix pencil $\mathbf{(}\tilde{\mathbf{E}}\mathbf{,}\tilde{\mathbf{A}}\mathbf{)}$ .

Figure 5

Output solutions of Example 4.

Example 5 In this example, we consider a RLC network descriptor model of electric power grid originating from [11]. It is called $\mathrm{M}\mathrm{8}\mathrm{O}\mathrm{P}\mathrm{I}$ system and can also be downloaded from [10]. It is an index-2 DAE system of dimension $n=4\text{,}182$ with 3 inputs and 3 outputs. The inputs are current signals injected into three buses given by:

$$u=\left[\begin{array}{c}sin({\omega }_1(t-t_1))U(t-t_1)\\ sin({\omega }_2(t-t_2))U(t-t_2)\\ sin({\omega }_3(t-t_3))U(t-t_3)\end{array}\right],$$

where $U(t)$ is the step function, ${\omega }_1=2\times \pi \times 2\text{,}000\mathrm{rad}/s$, ${\omega }_2=2\times \pi \times 750\mathrm{rad}/s$ and ${\omega }_3=2\times \pi \times 150\mathrm{rad}/s$. Values $t_1=5\mathrm{ms}$, $t_2=3\mathrm{ms}$ and $t_1=1\mathrm{ms}$ are time delays associated with the input signals. The desired solutions are voltages $e_1(t)$, $e_2(t)$ and $e_3(t)$ at the three buses. The sparsity of the matrix pencil of this DAE system is shown in Figure 6. We were able to decoupled the DAE system into $4\text{,}028$ differential equations and 154 algebraic equations using the procedure we discussed in Section 3.2.1, i.e. $n_{10}=4028$, $k_1=35$ and $k_0=119$. The sparsity of the matrix pencil of the decoupled system in the descriptor form is shown in Figure 7. We solved the differential part of the decoupled system using $\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{23}\mathrm{t}$ solver and we obtained good solutions as shown in Figure 8. The simulation took 39 seconds. We could not compare with the Matlab inbuilt solvers since they cannot solve index-2 systems. We validated our results with those presented in [11].

Figure 6

Sparsity of matrix pencil $\mathbf{(}\mathbf{E}\mathbf{,}\mathbf{A}\mathbf{)}$ .

Figure 7

Sparsity of matrix pencil $\mathbf{(}\tilde{\mathbf{E}}\mathbf{,}\tilde{\mathbf{A}}\mathbf{)}$ .

Figure 8

Output solutions of Example 5.

6 Conclusion

In conclusion, splitting DAEs using special projectors leads to decoupled systems which are easier to solve than DAE systems directly. We have tested the Split-DAE method on both large and small problems and proved to be a robust and efficient method. We can now use any numerical integration method to solve DAE systems. Hence this method can be used to simulate electromagnetic descriptor models.