US20090012770A1

US20090012770A1 - Circuit simulation

Info

Publication number: US20090012770A1
Application number: US12/060,556
Authority: US
Inventors: Oleg Wasynczuk; Juri V. Jatskevich
Original assignee: Pc Krause & Associates Inc
Current assignee: Pc Krause & Associates Inc
Priority date: 2001-01-11
Filing date: 2008-04-01
Publication date: 2009-01-08
Also published as: WO2002056145A2; AU2002243494A1; US7353157B2; WO2002056145A3; US20020183990A1

Abstract

A system, method, and apparatus select state variables for, build state equations of, and simulate time-domain operation of an electronic circuit. The circuit is modeled with three branch types (inductor, resistor, voltage source in series; capacitor, resistor, current source in parallel; and switch), including four pre-defined switch types (unidirectional unlatched, bidirectional unlatched, unidirectional latched, and bidirectional latched). Automated analyses determine efficient state variables based on the currently active circuit topology, and state equations are built and applied. Switching logic determines when switch states change, and state equations for the new topology are either drawn from a cache (if the topology has already been processed) or derived anew. The switch control signals may be combined into a single switching variable, defined as a function of the state output.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a continuation of U.S. patent application Ser. No. 10/043,981 filed Jan. 11, 2002, which claims the benefit of U.S. Provisional Patent Application 60/261,033, filed Jan. 11, 2001, both of which are incorporated herein by reference.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT

This invention was made with Government support under Contract F33615-99-C-2911 awarded by the U.S. Department of Defense. The United States Government has certain rights in the invention.

BACKGROUND

The present invention relates to the design, modeling, simulation, and emulation of electronic circuitry. More specifically, the present invention relates to numerical time-domain simulation of analog or digital electrical circuits using mathematical expressions. Present simulation systems suffer from limitations in the kinds and topologies of circuits to which they may be applied. The complexity of systems to be simulated is also limited in current systems by various inefficiencies in the simulation modeling process, including state selection and state equation building methods. There is thus a need for further contributions and improvements to circuit simulation technology.
Existing techniques for circuit simulation include two approaches. In one, the underlying numerical algorithms are based upon a state-space model of the system being simulated or analyzed. The system to be simulated is specified either in the form of a text file using a well-defined syntax, or graphically using boxes or icons to represent computational elements such as summers, multipliers, integrators, function generators, and/or transfer functions, for example. General purpose mathematical packages operate the models by assembling a state-space model of the system from the user-supplied description. The time-domain response is then calculated numerically using established integration algorithms.
In the other category of solutions, the system is described as an electrical circuit, the fundamental branches of which may include resistors, inductors, capacitors, voltage sources, and/or current sources, for example. Other circuit elements such as diodes, transistors, or electromechanical devices may also be defined as a variation or a combination (sub-circuit) of the fundamental branches. In this type of simulation system, the model developer describes the circuit in the form of a list of parameters of the fundamental or user-defined circuit elements and the layout of the circuit. Nodal or modified-nodal analysis techniques are then employed to simulate operation of the circuit. In that process, the differential equations associated with each inductive and capacitive branch are modeled in discrete form using difference equations to relate the branch current and voltage at a given instant of time to the branch current and/or voltage at one or more preceding instants of time. The difference equations for the overall circuit are then assembled automatically and solved using established methods to simulate the time-domain response.
It is of note that the second category of systems does not force the model developer to derive the state equations or block-diagram representation of the circuit. On the other hand, a state-space model of the overall circuit is never established or calculated by the program. Consequently, the numerous (and generally more efficient) techniques for analysis of state-space models cannot be applied in such systems.
More recently, a system has been developed that produces a state-space realization for the linear part of a system, the non-linear portion and variable parameters being represented as separate or external blocks. The overall system is then simulated using the appropriate solvers in the host mathematical package. In these systems, however, switches are modeled as resistors with very small or very large values. The resulting state model, therefore, may be artificially stiff and have larger dimension than necessary because of the states artificially introduced by the resistor-based switch model. In addition, the system does not incorporate mutual, non-linear, or variable inductances, capacitances, and resistances in its standard library blocks and components.
Some work on an automated state model generator and circuit simulator (ASMG) is reported in O. Wasynczuk and S. D. Sudoff, “Automated State Model Generation Algorithm for Power Circuits and Systems,” IEEE Transactions on Power Systems, Vol. 11, No. 9, November, 1996, pp. 1951-1956. In this work, circuits are specified and analyzed using fundamental branches as shown in FIG. 1. Each fundamental branch includes a switch, resistor r_br, inductor L_br, voltage source e_br, conductance g_br, capacitor P_br, and current source j_br. Each parameter can be fixed or time-varying, and ideal components can be modeled by setting the remaining parameters to zero. Given the parameters for each branch and the list of nodes that the branches connect, the ASMG generates a state-space model of the overall circuit. The state-space representation is calculated and solved using numerical methods. If a change in switching state occurs, the state model generator then recalculates the state-space model and establishes the appropriate initial conditions for the new topology. A disadvantage of this system is that it cannot be used to simulate circuits that include loops composed of voltage sources, capacitors, and/or resistors. This limitation dramatically hindered the ability of the ASMG to simulate high-frequency switching transients of power-electronic-based systems.
Subsequent improvements to the ASMG is the fundamental branches illustrated in FIG. 2 using the same notation for parameters as was used in FIG. 1. Using this form of modeling, more complicated circuit elements could be represented, including transistors, diodes, and thyristors. These improvements are discussed in J. Jatskevich, “A State Selection Algorithm for the Automated State Model Generator,” Ph.D. Thesis, Purdue University, 1999. Additional improvements are described herein.
As used herein, a “spanning tree” over a graph (comprising a set of branches and a set of nodes to which the branches are connected) is defined to be a subset of the set of branches such that at least one branch in the subset is connected to each node from the set of nodes, yet no loop can be formed from the branches in the subset.
Also, as used herein, a “topology change event” occurs when the control signal for one or more switching elements causes the switching state of that element to change. The nodes to which that element is connected will then become (or cease to be) galvanically connected to another portion of the circuit. In most cases, this change affects the state equations for the overall circuit.

SUMMARY

It is an object of the present invention to provide an improved system, method, and apparatus for simulating electrical and electronic circuits. Another object is to provide a system, method, and apparatus to more efficiently simulate the operation of a wider variety of electrical and electronic circuits than has previously been possible.
These objects and others are achieved by various forms of the present invention. One form of the present invention is a method, including creating one or more data structures sufficient to model an electronic circuit as a collection of n (at least two) elements. These comprise zero or more LRV elements, zero or more CRI elements, and zero or more switching elements. The LRV elements each have at least one of (a) a non-zero inductance parameter L_br, (b) a non-zero resistance parameter r_br, and (c) a non-zero voltage source parameter e_br, but neither a non-zero capacitance parameter, nor a non-zero current source parameter, nor a switch parameter. The CRI elements each have at least one of (a) a non-zero capacitance parameter C_br, (b) a non-zero resistance parameter r_br, or (c) a non-zero current source parameter j_br, but neither a non-zero inductance parameter, nor a non-zero voltage source parameter, nor a switch parameter. The switching elements each have a switch state and neither a non-zero inductance parameter, a non-zero capacitance parameter, a non-zero resistance parameter, a non-zero voltage source parameter, nor a non-zero current source parameter. A first set of state equations is automatically generated from the one or more data structures, and operation of the electronic circuit is simulated by application of said first set of state equations. In this method, the collection comprises either (1) an LRV element for which at least two of L_br, r_br, or e_brare non-zero, or (2) a CRI element for which at least two of C_br, r_br, or j_brare non-zero.
In variations of this form of the invention, the simulating step includes producing state output data. In some embodiments, some or all of the parameters in the first set of state equations change over (simulation) time as a function of the state output data. In some embodiments, some or all of the parameters change over (simulation) time due to a time-varying parameter of at least one element in the collection.
In other variations of this form of the invention, a second set of state equations is generated from the one or more data structures upon the occurrence of a first topology change event. In some such embodiments, the generating step simply involves modifying only the subset of said first set of state equations that depends on the one or more switching elements that have changed. In other such embodiments, each unique vector of switch states represents a topology of the overall circuit, and the method also includes (1) storing the first set of state equations in a cache; (2) after a second topology change event, determining whether a set of state equations in the cache represents the new topology; and (3a) if the determining step is answered in the affirmative, using the set of state equations that represents the new topology to simulate operation of the circuit after the second topology change event; or (3b) if the determining step is answered in the negative, building a third set of state equations that represents the new topology, and using the third set of state equations to simulate operation of the circuit after the second topology change event.
In other such embodiments, the method also includes (1) storing said second set of state equations in a cache (2) after a third topology change event, deciding whether a set of state equations in the cache represents the new topology; and (3a) if the deciding step is concluded in the affirmative, using the set of state equations from the cache that represents the new topology to simulate operation of the circuit after the third topology change event; or (3b) if the deciding step is concluded in the negative, building a new set of state equations that represents the new topology, and using the new set of state equations to simulate operation of the circuit after the third topology change event.
Another form of the present invention is a method including creating one or more data structures that together store characteristics of some active branches B^activethat make up a graph of nodes and branches that form a circuit, wherein B^activeconsists of (i) a (possibly empty) set B^Lof inductive branches, each having a non-zero inductive component but neither a capacitive component nor a variable switch state; (ii) a (possibly empty) set B^Cof zero or more capacitive branches, each having a non-zero capacitive component but neither an inductive component nor a variable switch state; and (iii) a (possibly empty) set B^Aof additional branches, each having neither an inductive component nor a capacitive component. B^activeis partitioned into a first branch set B_tree ^activeand a second branch set B_link ^active, where the branches in B_tree ^activeform a spanning tree over B^active, giving priority in said partitioning to branches not in B^Lover branches in B^L. B_link ^activeis sub-partitioned into a third branch set B_link ^L, and a fourth branch set B_link ^CA, where B_link ^L=B_link ^active∩B^L. A fifth branch set B^CAis identified as the union of (i) B_link ^CA, (ii) B^C∩B_tree ^active, and (iii) those branches in B_tree ^activethat form a closed graph when combined with B_link ^CA. B^CAis partitioned into a sixth branch set {tilde over (B)}_tree ^CAand a seventh branch set {tilde over (B)}_link ^CA, where the branches in {tilde over (B)}_tree ^CAform a spanning tree over B^CA, giving priority in said partitioning to branches in B^Cover branches not in B^C. An eighth branch set B_tree ^C={tilde over (B)}_tree ^CA∩B^Cis identified. A set of state variables is selected, comprising (a) for each branch of B_link ^L, either the inductor current or inductor flux, and (b) for each branch of {tilde over (B)}_tree ^C, either the capacitor voltage or capacitor charge. A plurality of states of the circuit are simulated using the set of state variables.
In a variation on this form of the invention, the partitioning steps each comprise an application of a weighted spanning tree algorithm, such as, for some positive numbers w_Land w_C, (a) for the partitioning of B^active, a minimum spanning tree algorithm is used with weight function
$ω_{L} (b_{j}) = {\begin{matrix} w_{L} if branch b_{j} \in B^{L} \\ 0 otherwise; \end{matrix}$
and (b) for the partitioning of B^CA, a maximum spanning tree algorithm is used with weight function
$ω_{C} (b_{j}) = {\begin{matrix} w_{C} if branch b_{j} \in B^{C} \\ 0 otherwise . \end{matrix}$
Another form of the invention is a system, including a processor and a computer-readable medium in communication with the processor, where the medium contains programming instructions executable by the processor to:
(1) build state equations for a first topology of an electronic circuit that has at least two switching elements, where each switching element has a switching state;
(2) solve the state equations at some time t_ito provide a state output vector, in which at least two elements control the switching states of the switching elements;
(3) calculate the value of a switching variable as a function of the state output vector, wherein the value reflects whether the switching state of at least one of the switching elements is changing; and (4) if the value of the switching variable at time t_iindicates that at least one of the switching elements is changing, determine a second topology of the electronic circuit for time t_i ⁺ and obtain state equations for the second topology.
In a variation of this form of the invention, the programming instructions include a state equation building module, a solver module for ordinary differential equations, and a switching logic module. The building is performed by the state equation building module, the solving and calculating are performed by the solver module; and the determining is performed by the switching logic module. In some embodiments, the obtaining is performed by said switching logic module, while in others the obtaining is performed by said state equation building module.
In still other embodiments, at some time t_j, at least two switching elements are each either rising-sensitive or falling-sensitive switches. Rising-sensitive switches change switching state if and only if a controlling element of the state vector has passed from a negative value to a non-negative value. Conversely, falling-sensitive switches change switching state if and only if a controlling element of the state vector has passed from a positive value to a non-positive value. In these embodiments, the function is the arithmetic maximum of (1) a maximum of all elements of the state vector that control rising-sensitive switches, and (2) the negative of the minimum of all controlling elements of the state vector that control falling-sensitive switches.
A still further form of the invention is a system for simulating electronic circuits having a processor and a computer-readable medium in communication with said processor, where the medium contains programming instructions executable by the processor to read element parameters and node connection information from a data stream. The stream includes at least one switch type specification selected from the group consisting of: a unidirectional, unlatched switch; a bidirectional, unlatched switch; a unidirectional, latched switch; and a bidirectional, latched switch. The instructions are further executable by the processor automatically to calculate state equations for the circuit given the states of switches specified by the at least one switch type specification.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram of a generic branch for modeling a circuit using a prior art system.

FIG. 2 is a diagram of two branches used to model circuits in another prior art system.

FIG. 3 is a block diagram of the overall process for simulating a circuit according to the present invention.

FIGS. 4A-4C are branches used to model circuit components using the present invention.

FIG. 5 is a block diagram of the run-time computational routines in one embodiment of the invention.

FIG. 6 is a block diagram detailing the inductive link current calculator for use with the routines shown in FIG. 5.

FIG. 7 is a block diagram detailing the capacitive tree voltage calculator for use with the routines shown in FIG. 5.

FIG. 8 is a block diagram detailing the resistive network algebraic equation calculator for use with the routines shown in FIG. 5.

FIG. 9 is a block diagram detailing the inductive network state/output equation calculator for use with the routines shown in FIG. 5.

FIG. 10 is a block diagram detailing the capacitive network state/output equation calculator for use with the routines shown in FIG. 5.

FIG. 11 is a schematic diagram of a switch element.

FIG. 12 is a block diagram of the interaction between subnetworks as analyzed for use with the present invention.

DESCRIPTION OF THE ILLUSTRATED EMBODIMENTS

For the purpose of promoting an understanding of the principles of the present invention, reference will now be made to the embodiment illustrated in the drawings and specific language will be used to describe the same. It will, nevertheless, be understood that no limitation of the scope of the invention is thereby intended; any alterations and further modifications of the described or illustrated embodiments, and any further applications of the principles of the invention as illustrated therein are contemplated as would normally occur to one skilled in the art to which the invention relates.
Generally, the system illustrated in FIGS. 3-10 and described herein simulates operation of a circuit in the time domain by collecting component parameters (variable and/or constant) and the overall circuit topology, then establishing a minimal state space for each active topology (as they are encountered), building state equations for that state space, and solving those equations for relevant steps in time. As discussed herein, the above process can be very much streamlined using the additional knowledge and techniques provided in the present invention.
The overall simulation process will now be discussed with reference to FIG. 3. A description 21 of the circuit, including an identification of constant and variable parameters, initial conditions, and the like are fed into state model generator 31, which provides state model 23 to solver 33. Solver 33 generates the simulation out put 25 for data consumers (such as log files, graphical visualization tools, and the like) as are known in the art. Solver 33 also provides continuous state information 27 to switching logic 35, which determines the state of one or more switches in the circuit should be changed. The result of this analysis, switching state 29, is provided to state model generator 31. If topological changes to the circuit are indicated in switching state 29, state model generator 31 updates the state model and passes that updated model to solver 33 for continued simulation.

Finite Electrical Networks

For the purposes of automated modeling, it is assumed that electrical networks can be composed of inductive, capacitive, and switch branch models depicted in FIGS. 4( a), (b), and (c), respectively. Using this approach, a wide variety of electrical circuits with different topologies can be modeled by appropriately setting the branch parameters. Such a modeling approach also assumes that only a finite number of branches is allowed for representation of any particular circuit, and that all branches have lumped parameters. Electrical networks satisfying this assumption can be modeled using a finite-dimensional state variable approach. Thereafter, it is possible to derive a finite-dimensional system of ordinary differential equations (ODEs) and algebraic equations (AEs) that would portray the dynamic behavior of currents and voltages for any such circuit. It is also important to note that only “practical” or “reasonable” circuits are considered, which means finite energy, finite current, and finite voltage electrical networks. For such networks, all energy storage elements such as capacitors and inductors are allowed to store only a finite amount of energy at any instant of time and that energy cannot change forms instantaneously. Also, no branches are allowed to carry infinite current through their components or have infinite voltage between nodes throughout the existence of the network including the instant of time that the network is switching between topological stages. These restrictions ensure a proper commutation (or proper transition) of the network from one topology to another. Therefore, a class of electrical networks that permits only proper commutation and posses a finite dimensional system of differential algebraic equations (DAEs) can be defined as a class of finite networks. Such a class of finite electrical networks is denoted as N_q ^α+β where α and β are the number of state variables in the systems of ODEs for the inductive and capacitive networks, respectively, and q is the number of branches used to represent the layout of the network. Also, α and β are often called the network complexity. Thus, hereafter only proper and finite electrical networks from the class N_a ^α+β are considered.

Representation of Electrical Networks

This section will define certain terms with respect to electrical networks from the class N_a ^α+β. In particular, a network N can be defined by its topological configuration, which is best described by the associated graph denoted as G, and a set containing branch parameters denoted as P. Also, a particular switch branch may be identified as active or inactive depending on its state. Therefore, in order to specify N completely, a topological vector s, which would contain information regarding whether each branch is active or inactive, should be added. A topological state vector s would have ones in places corresponding to all currently active branches, and zeros for the remaining branches which were identified as inactive for the current topology. Thus a network of a general kind is an object from the class N_a ^α+β and is defined as a triple
N=(G,P,s) (2.1)
It is sometimes convenient not having to deal with inactive branches and assume that all of the network branches are currently active. Such an assumption would significantly simplify notation without a significant loss of generality. Often, a network with all of its branches assumed to be active may be referred to augmented or having augmented topology. In these cases, the vector s would not carry any additional information and, therefore, can be omitted from (2.1).
For the associated graph G to be a circuit, it is required that any single branch of the overall network N must be included in some loop with other branches so as to ensure the existence of a closed path for the current flow. Graphs with such a circuit property are referred to as being closed. For consistency, a reference direction (from − to +) is assigned to each branch with respect to the two nodes at its endpoints. Also, a network N, in the most general case, may be composed of several galvanically independent circuits. Therefore, a directed multi-graph denoted by G_d ^gcan be associated with network N. Such a multi-graph may consist of g closed subgraphs. That is
G_d ^g={G_d ¹, G_d ², . . . G_d ^k, . . . G_d ^g} (2.2)
where each k-th closed subgraph G_d ^kcorresponds to its circuit. The subscript d indicates that the graph is directed. Even thought all such subgraphs are electrically disjointed they may be coupled through mutual inductances and/or capacitances between branches. On the other hand, if g=1, the associated network becomes a single circuit. Sometimes, for the purposes of analysis, it is also necessary to consider the associated undirected graph or subgraph, which can be denoted at G^gand G^k, and obtained by omitting the information regarding the reference direction of each branch. In a closed undirected graph, any branch must be a member of some cycle with more than one branch.
The directed multi-graph G_d ^gconsists of a total of q branches from the branch set
B={b₁, b₂, . . . b_j, . . . , b_q} (2.3)
which, in turn, are connected in some arrangement via p nodes from the node set
N={n₁, n₂, . . . , n_i, . . . n_p} (2.4)
Branches and nodes appear in (2.3) and (2.4) in a definite order which is given by their respective subscript indices. Using such a representation, it is possible to retrieve a particular branch or a node by referring to the respective ordered set, B or N, with an appropriate index-pointer.
Therefore, the graph is defined in terms of the branch and node sets as G_d ^g=(N,B). In general, the network N has a changing topology. At one instant of time, some branches may be inactive and some nodes may be unused. At another instant of time, other branches may switch on and off including different nodes, thus defining a new topology. Whenever it is needed, the number of currently active branches may be denoted as q′ and the number of currently used nodes as p′; whereupon the number of inactive branches and nodes may be expressed as q-q′, and p-p′, respectively. Often, when it is clear what the node and branch sets are, only their dimensions may be specified instead of actual sets. Also, a spanning tree G_tree ^kcan be associated with each k-th closed subgraph G^k, and a corresponding forest of such trees G_tree ^g=(N,B_y) can be associated with the global graph G^g=(N,B). Note that the set of forest branches B_yis the subset of the global branch set, that is B_y⊂B.
There are several ways to represent a graph on a computer. Some methods may take advantage of the sparsity of the interconnection matrix and are, therefore more efficient than others in terms of the memory required to store a particular graph. A less efficient but algebraically convenient method is to represent a graph in matrix form. In particular, a node incidence matrix A_ffor the multi-graph has p rows and q columns (one row for each node and one column for each branch, all ordered). Even though this matrix never has full rank, it can be referred to as a full node incidence matrix meaning that it is not reduced to the size of its rank. This matrix is conveniently formed from positive and negative node incidence matrices as
A _f =A ₊ −A ₋ (2.5)
where the entries of A₊ are A₊(i, j)=1 if the j-th branch includes the i-th positive node, and zero otherwise. Similarly, the only non-zero entries of A₋ are A₋(i, j)=1 if the j-th branch includes the i-th negative node. For the networks with changing topology, this matrix can be updated for each topology such that for each currently inactive branch b_jthe corresponding j-th column of A_fis replaced with zeros. The transpose of A_fis also known as the adjacency matrix of a graph. For large graphs, these matrices are quite sparse and, therefore, may be stored using techniques optimized for sparse matrices.
It is understood in the art that each connected subgraph of the network N contributes exactly one zero-row to A_RREF. Thus, the rank of A_fis always the number of active nodes less the number of connected subgraphs, which can be written as
rank(A _f)=rank(A _RREF)=p′−g=r (2.8)
It can be observed that the total number of loops in the undirected graph G^g=(N,B) associated with N may be computed as the number of active branches less the rank of A_f, which may be written as
s=q′−rank(A _f)=q′−r (2.9)
Using the branch models in FIG. 4, the parameters of the circuit can be conveniently arranged into parameter matrices denoted as R_br, G_br, L_br, and C_br. In one embodiment, these square matrices would have dimensions equal to the total number of branches in the global network to be modeled, and would contain the parameters of each branch in the corresponding diagonal entry of each matrix. If some branches are coupled through mutual inductances or mutual capacitances, then those mutual parameters are represented in the off-diagonal entries of L_bror C_brcorresponding to these branches. In the same way, mutual resistances or mutual conductances can be modeled between branches. Even though such quantities might not have physical meaning similar to mutual inductances or mutual capacitances, they can be employed for simulation purposes. If some branch parameters are not present in the circuit, the corresponding entries in the parameter matrices are filled with zeros. Similar to parameter matrices, the external independent voltage and current sources are also represented as vectors denoted as e_brand j_br. These vectors would have voltages and currents in those entries corresponding to the branches with respective sources and zeros elsewhere.
In some circuits to be simulated, the network parameters will not depend on the currents and voltages applied to the branches. In other circuits the parameters will vary with time, so the parameter matrices would also become matrix functions of time R_br(t), G_br(t), L_br(t), and C_br(t). For simplicity of notation herein, such time dependence will not be written explicitly. Also, since inductors and capacitors are energy storage components, in addition to their time dependence, their total derivatives with respect to time must also be known. Thus, the derivatives of time-varying inductances and capacitances are additional inputs into the model of the energy transformation and exchange processes in the circuit.
The network parameter matrices and the vectors of voltage and current sources can be grouped to form a parameter set for the global network, which can be defined as
$\begin{matrix} P = {R_{br}, L_{br}, \frac{\partial L_{br}}{\partial t}, e_{br}, G_{br}, C_{br}, \frac{\partial C_{br}}{\partial t}, j_{br}} & (2.10) \end{matrix}$
It is also convenient to define the following subsets of P
$\begin{matrix} P_{L} = {R_{br}, L_{br}, \frac{\partial L_{br}}{\partial t}, e_{br}} & (2.11) \\ P_{C} = {G_{br}, C_{br}, \frac{\partial C_{br}}{\partial t}, j_{br}} & (2.12) \\ P_{A} = {R_{br}, G_{br}, e_{br}, j_{br}} Clearly, & (2.13) \\ P = P_{L} ⋃ P_{C} ⋃ P_{A} & (2.14) \end{matrix}$
Subdivision of the parameter set P corresponding to the global network N into subsets (2.11), (2.12), and (2.13) has the following goal. First, consider a grid of branches of either of the two types such that all parameters can be placed in one set P_A. A network of this kind does not have energy storage element such as inductors or capacitors. Modeling such a network would require solving a system of AEs relating branch voltages and currents to external sources through the network parameters and topology. From the nature of the equation that needs to be solved in order to obtain all branch currents and voltages, such networks are referred to as being algebraic and denoted as N_A. On the other hand, networks whose parameters are placed in the set P_Lhave inductive elements which may store energy in the magnetic field. Therefore, a network of this type possesses a system of ODEs, or more precisely a state equation, whose natural state variable may be inductor currents or flux linkages. Similarly, networks whose parameters are placed in PC have capacitors that may store energy in an electric field, and therefore, posses ODEs whose natural state variables may be capacitor voltages or charges. These two networks are referred to as inductive and capacitive, respectively.
Thus, three types of elementary electrical networks, namely inductive, capacitive, and algebraic, denoted as N_L, N_C, and N_A, respectively, are considered. There, only the first two networks are allowed to have energy storage circuit components, and therefore, only these two networks can have state variables. In the most general case, an electrical system may be composed of all types of branches considered above. Therefore, a corresponding network N can be represented as an interconnection of N_L, N_C, and N_A, which is symbolically expressed in terms of union as
N=N_L∪N_C∪N_A (2.15)
Depending on a particular circuit to be modeled, some networks may or may not be present in (2.15).

Topological Canonical Form

Models of electrical networks must obey all laws of circuit theory, and in particular, Kirchhoff's current and voltage laws, KCL and KVL. Automated modeling requires formulation of equations describing the network using KCL and KVL that are written algorithmically for the entire circuit or its sections. Moreover, KCL and KVL may be written for a network of the most general kind based only on its topological configuration and regardless of the actual branch elements of their volt-ampere characteristics. Therefore, the starting point in the analysis of graph G_d ^g=(N,B) is the associated full node incidence matrix A_f.
It is necessary to express KCL and KVL for the network whose topology is given by G_d ^g=(N,B) and stored in A_f. For simplicity of notation, all branches in this section are assumed to be active. The procedure continues as follows. First, instead of G_d ^g=(N,B), the corresponding undirected version G^g=(N,B) is considered. Then, suppose that from G^g=(N,B), it is possible to find a forest of spanning trees G_trees ^g=(N,B_y) that spans the entire node set N using the set B_y, which is some subset of B of minimal size over all branch subsets. In general, G_trees ^g=(N,B_y) is far from being unique for a given G^g=(N,B). In fact, it can be shown that a graph whose full node incidence matrix is A_fhas det(A_fA_f ^T) total number of spanning trees. Moreover, the partitioning of a multi-graph into spanning trees G_trees ^g=(N,B_y) and remaining link-branches is a nondeterministic topological problem. One way to simplify this task is to associate appropriate weights with the branches and convert the problem into finding a spanning tree with the minimized/maximized sum of such weights. This approach is known in network optimization as the minimum/maximum spanning tree problem. There are several well-known minimum/maximum spanning tree algorithms which roughly yield the order of complexity O(q log p) and better depending on the data structure used to represent a graph. Also, the branch weights should be relatively simple so as to promote good performance as well as to be able to prove certain useful properties. Thus, the weights can be assigned to network branches based on their respective parameters (2.10). This method of obtaining a spanning tree and a set of links, with some desired property, will be utilized in the present, exemplary embodiment for the purpose of automated modeling.
For the purpose of this section, any proper forest of spanning trees would suffice. Having performed such graph partitioning, the set of tree-branches B_yis identified. Thereafter, the set of remaining link-branches B_xcan be determined as
B _x =B−B _y (2.16)
Since the subsets B_yand B_xare identified based on G_trees ^g=(N,B_y), it is now also possible to re-sort the complete branch set B such that all tree-branches appear first from the left and all the remaining link-branches on the right. That is, the branch set B can be re-ordered as
B={B_y,B_x} (2.17)
The new branch order defined in (2.17) can be applied to the columns of A_f. In this case, the new branch order can be related to the original branch order through a permutation matrix which is denoted as T_p. Multiplying A_ffrom the right by T_presults in a matrix whose columns are ordered to correspond to the branch set B in (2.17). That is, the permutation matrix T_pshould sort the columns of A_fsuch that
A_fT_p=[A_trees,A_links]=Â_f (2.18)
This permutation matrix can be assembled from an identity matrix by sorting its columns at the same time as the branches in (2.17). Note that the multiplication from the right by T_p ^Tperforms the reverse column permutation and restores the original column order. That is
A_f=Â_fT_p ^T (2.19)
Here, as well as further on, the hat sign above the matrix denotes that the corresponding matrix or vector quantity is referred to a branch order different that the original order.
After the full node incidence matrix is expressed in the form (2.18), its algebraic properties can be utilized. In particular, since the columns of A_treescorrespond only to branches that form a free-of-cycles forest, this matrix must have a full rank. On the other hand, A_linkscontains only columns corresponding to the link-branches, and therefore, none of its columns add to the rank of Ã_f. Therefore, the rank of Ã_f, as well as the rank of A_f, is determined by the columns of A_trees. Therefore, these columns can be chosen to be the basis columns of Ã_for A_f. With respect to (2.18), the basis columns appear first from left-to-right. The RREF of (2.18), therefore, will have the identity matrix of the size of its rank as the upper-left block, and its re-ordered structure becomes
$\begin{matrix} RREF ({\hat{A}}_{f}) = [\begin{matrix} I_{rxr} {\hat{A}}_{rxs} \\ 0_{gxr} 0_{gxs} \end{matrix}] & (2.20) \end{matrix}$
The zero rows on the bottom of (2.20) do not include any useful topological or algebraic information about the network and therefore can be deleted. We define the topological canonical form (TCF) of the node incidence matrix to be the result of the reduction of (2.20):
A _TCF=TCF(A _f)=└I _rxr Â _rxs┘ (2.21)
The TCF is one of the key concepts used herein to describe networks. From A_TCF, the reduced incidence matrix Â_aand the so called basic loop matrix {circumflex over (B)}_bare obtained as
Â_a=A_TCF (2.22)
{circumflex over (B)} _b =└−Â _sxr ^T I _sxs┘ (2.23)
Matrices (2.22) and (2.23) were assembled with respect to the branch order (2.17) and can be easily recovered for the original branch order as
A_a=Â_aT_p ^T (2.24)
B_b={circumflex over (B)}_bT_p ^T (2.25)
Denoting vectors of branch currents and voltages as i_brand v_br, respectively, KCL and KVL written for the whole network have the usual form
A_ai_br=0 (2.26)
B_bv_br=0 (2.27)
Also, instead of using recovered matrices (2.24)-(2.25), it is possible to re-sort i_brand v_brinto î_brand {circumflex over (v)}_brfor the new branch order (2.17) using the same permutation matrix T_p, and then use these with (2.22)-(2.23) to form KCL and KVL similar to (2.26)-(2.27). That is
Â_ai_brT_p=Â_aî_br=0 (2.28)
Â_av_brT_p={circumflex over (B)}_b{circumflex over (v)}_br=0 (2.29)
Furthermore, based on the structure of TCF, an interesting property of matrices (2.22) and (2.23) can be utilized. This property relates voltages and currents of the branches corresponding to subsets B_yand B_xin (2.17). In particular, based on the branch order (2.17) it is possible to define the vectors of tree and link currents and voltages as
î_br=[i_y,i_x] (2.30)
{circumflex over (v)}_br=└v_y, v_x┘ (2.31)
Then, based on (2.22), (2.23), and (2.30), the vector of re-sorted branch currents can again be expressed in terms of the link currents as
î _br=({circumflex over (B)} _b)^T i _x (2.32)
Similarly, the vector of re-sorted branch voltages can be expressed in terms of the tree voltages as
{circumflex over (v)} _br(Â _a)^T v _y (2.33)
The relations (2.28)-(2.33) play an important role in formulating the governing DAEs describing electrical networks.
Networks with State Variables
In this section, a class of finite electrical networks N_q ^α+β is considered. In particular, the two types of networks that can be modeled using a state variable approach are inductive and capacitive. Equipped with techniques based on a topological search for an appropriate forest of spanning trees, the conditions under which a corresponding state equation can be assembled are also considered.

Inductive Network State Equation Formulation

An inductive network can be built using branches of the type depicted in FIG. 4( a). Assuming that all switches are active, a network of this type is given as (2.80)
N _L=(G,P _L) (2.80)
where the parameter set P_Lis defined in (2.11).
For a given topology of N_L, the branches must be re-sorted into subsets as in (2.17). The subset B_ywould collect branches that form spanning trees for all subgraphs. Each such tree is free of cycles, and covers all active nodes in its subgraph. The second category, denoted as subset B_x, takes the link-branches. These branches are the links in a sense that addition of any of them to the spanning tree would result in a cycle. These branches can carry state variables—independent currents—and since their number is minimal for each spanning tree, they form a minimal set of states for the N_L. Thus, in order to partition G^g=(N,B) and obtain the required branch sets, the following weight function w_L(b) is defined as
$\begin{matrix} For j = 1 \dots q, w_{L} (b_{j}) = {\begin{matrix} 0, if L_{br} (j, j) = 0 \\ 1, if L_{br} (j, j) \neq 0 \end{matrix} & (2.81) \end{matrix}$
Applying the minimum spanning tree algorithm (MinSTA) to G^g=(N,B) with weight function w_L(b) a minimum spanning forest G_trees ^g=(N,B_y) is obtained as before. Thereafter, the set of link branches B_xcan be determined as in (2.16).
An important observation is that the network N_Lhas no non-inductive loops if and only if there is a set B_xin which the number of branches equals the sum of their weights
$\begin{matrix} \sum_{j \in B_{x}} w_{L} (b_{j}) = size (B_{x}) = s & (2.82) \end{matrix}$
Condition (2.82) in necessary and sufficient. The set B_xneed not be unique, but any such set B_x, satisfying (2.82) is equivalent in a topological sense. That is, an arbitrary set of branches (that can be larger than B_x) may be chosen to represent state variables, such as independent currents, in N_Lif and only if it contains a set B_xsatisfying topological condition (2.82). If however, (2.82) cannot be met, it follows that the given network is more than just a single inductive network. An algorithm for handling more than one network will be presented in later sections.
Applying the MinSTA with (2.81), assuming (2.82) holds, the order in which all branches are grouped and sorted according to (2.17) is obtained. This final branch order is related to the original branch order through a permutation matrix, which in this case is denoted as T_L. Using this permutation matrix and the TCF, matrices A_aand B_bare found as usual. Thereafter, a vector of state variables for N_Lmay be chosen to be a vector of independent currents, such that
i_br=B_b ^Ti_x (2.83)
The corresponding state equation is obtained using the dimensionality reduction procedure discussed in Wasynczuk and Sudhoff. The procedure may be as follows. The voltage equation written for the network is multiplied from the left by the corresponding matrix of KVL and all vectors of branch currents are replaced with (2.83). The result is
$\begin{matrix} \frac{\partial i_{x}}{\partial t} = - L_{x}^{- 1} (R_{x} + \frac{\partial L_{x}}{\partial t}) i_{x} - L_{x}^{- 1} B_{b} e_{br} & (2.84) \end{matrix}$

Capacitive Network State Equation Formulation

A capacitive network with augmented topology can be defined using the corresponding parameter set and an associated graph as
N _C=(G,P _C) (2.85)
where the set P_Cis defined in (2.12).
Also, in the case of N_C, the forest B_yof the associated G_trees ^g=(N,B_y) can represent state variables the capacitor voltages. Thus, a slightly different topological approach might be used. In particular, a MaxSTA with a different weight function can be applied in order to find maximum spanning trees. This time, the weight function w_C(b_j) is defined such that
$\begin{matrix} for i = 1 \dots p w_{C} (b_{j}) = {\begin{matrix} 0, if C_{br} (j, j) = 0 \\ 1, if C_{br} (j, j) \neq 0 \end{matrix} & (2.86) \end{matrix}$
In this case, the network N_Cshould not have non-capacitive tree-branches. Such a condition is satisfied if and only if there is a branch set B_yfor which the following is true
$\begin{matrix} \sum_{b_{i} \in B_{y}} w_{c} (b_{i}) = size (B_{y}) = r & (2.87) \end{matrix}$
Condition (2.87) is also necessary and sufficient, and therefore, the discussion thereof applies to this case as well. Again, the set B_yneed not be unique, but all such sets for which (2.87) holds are equivalent in a topological sense. If the condition (2.87) cannot be met, it can be shown that the corresponding network is not just N_Cbut a union of the form N_C∪N_A. However, for the purpose of this section, a single capacitive network is considered.
Thus, applying the MaxSTA with weight function w_C(b) and condition (2.87), a similar permutation matrix T_C, can be formed such that multiplying A_ffrom the right by T_Cthe desired branch order (2.17) for the N_Ccan be obtained. Further, matrices A_a, and B_bare formed exactly as before.
The natural state variables for N_Care the capacitor voltages. Therefore, a vector of state variables v_ycan be chosen to be a vector of independent capacitor voltages such that
v_br=A_a ^Tv_y (2.88)
The state equation is very similar to (2.84) and is also obtained using the dimensionality reduction procedure applied to the current equation written in matrix form for N_C. The result is
$\begin{matrix} \frac{\partial v_{y}}{\partial t} = - C_{y}^{- 1} (G_{y} + \frac{\partial C_{y}}{\partial t}) v_{y} + C_{y}^{- 1} A_{a} j_{br} & (2.89) \end{matrix}$
where notations are analogous to those used in Wasynczuk and Sudhoff, supra. Specifically
$\begin{matrix} G_{y} = A_{a} G_{br} A_{a}^{T} & (2.90) \\ C_{y} = A_{a} C_{br} A_{a}^{T} & (2.91) \\ \frac{\partial C_{y}}{\partial t} = A_{a} (\frac{\partial C_{br}}{\partial t}) A_{a}^{T} & (2.92) \end{matrix}$

Transformation of Variables

An important observation regarding the selection of state variables for N_Las well as for N_Cis that, once the set of states in (2.84) and (2.92) are obtained based on topological conditions specified previously, none of the actual states are required to be associated with any particular branch. That is, a topologically proper set of state variables in N_Land N_Ccan be transformed into an equivalent set by any non-singular coordinate transformation. Such a transformation of states is nothing more than a change of variables in the state space. To show this, it is useful to introduce time-varying coordinate transformation defined by two non-singular square matrices K_Land K_Cof dimensions corresponding to the number of states in N_Land N_C, respectively. Multiplying (2.26) and (2.27) from the left by K_C ^Tand K_L ^Trespectively, KCL and KVL still hold:
K_C ^TA_ai_br=0 (2.99)
K_L ^TB_bv_br=0 (2.100)
Thereafter, the new state variables for N_Land N_Care related to i_xand v_yas
i_x=K_Lx_L (2.101)
v_y=K_Cx_C (2.102)
And the corresponding branch currents in N_Land the branch voltages in N_Care expressed as
i_br=B_b ^TK_Lx_L={tilde over (B)}_b ^Tx_L (2.103)
v_br=A_a ^TK_Cx_C=Ã_a ^Tx_C (2.104)
Since the change of variables for the two networks is very similar, only the corresponding state equations for N_Lwill be written explicitly. Proceeding as follows
$\begin{matrix} \begin{matrix} \frac{\partial x_{L}}{\partial t} = - K_{L}^{- 1} L_{x}^{- 1} (R_{x} + \frac{\partial L_{x}}{\partial t}) K_{L} x_{L} - K_{L}^{- 1} \frac{\partial}{\partial t} (K_{L}) x_{L} - \\ K_{L}^{- 1} L_{x}^{- 1} B_{b} e_{br} \\ = - (K_{L}^{- 1} L_{x}^{- 1} R_{x} K_{L} + K_{L}^{- 1} L_{x}^{- 1} \frac{\partial}{\partial t} (L_{x}) K_{L} + K_{L}^{- 1} \frac{\partial}{\partial t} (K_{L})) \\ x_{L} - K_{L}^{- 1} L_{x}^{- 1} B_{b} e_{br} \end{matrix} & (2.105) \end{matrix}$
Choosing K_Lto be a constant matrix, the result of transformation is
$\begin{matrix} \begin{matrix} \frac{\partial x_{L}}{\partial t} = - K_{L}^{- 1} L_{x}^{- 1} (R_{x} + \frac{\partial L_{x}}{\partial t}) K_{L} x_{L} - K_{L}^{- 1} L_{x}^{- 1} B_{b} e_{br} \\ = - {\tilde{L}}_{x}^{- 1} ({\tilde{R}}_{x} + \frac{\partial {\tilde{L}}_{x}}{\partial t}) x_{L} - {\tilde{L}}_{x}^{- 1} {\tilde{B}}_{b} e_{br} \end{matrix} & (2.106) \end{matrix}$
where the reduced quantities with the subscript “x” and tilde above are defined using {tilde over (B)}_bfrom (2.106) instead of just B_b. Since (2.84) and (2.106) are related through the similarity transformation via constant matrix K_L, the eigenvalues of the resulting dynamic matrices in both state equations are the same. Also, since it is possible to choose K_Lto be a permutation matrix, it follows that the branch ordering does not change the eigenvalues of the system.
It is also interesting to observe what happens to (2.105) if the time-varying coordinate transformation is chosen to be
K_L=L_x ⁻¹ (2.107)
In this case, the new state variables become fluxes
λ_x=K_L ⁻¹i_x=L_xi_x (2.108)
Then, with the following equality
$\begin{matrix} \frac{\partial}{\partial t} (K_{L}^{- 1}) K_{L} + K_{L}^{- 1} \frac{\partial}{\partial t} (K_{L}) = 0 & (2.109) \end{matrix}$
equation (2.105) reduces to state equation
$\begin{matrix} \frac{\partial λ_{x}}{\partial t} = - R_{x} L_{x}^{- 1} λ_{x} - B_{b} e_{br} & (2.110) \end{matrix}$
A similar change of variables can be performed on N_C. The equations can be simplified even further by considering only constant capacitances, wherein a time invariant coordinate transformation K_Ccan be used. Therefore, (2.105) may be rewritten as
$\begin{matrix} \frac{\partial x_{C}}{\partial t} = - K_{C}^{- 1} C_{y}^{- 1} G_{y} K_{C} x_{C} + K_{C}^{- 1} C_{y}^{- 1} A_{a} j_{br} = - {\tilde{C}}_{y}^{- 1} {\tilde{G}}_{y} + {\tilde{C}}_{y}^{- 1} {\tilde{A}}_{a} j_{br} & (2.111) \end{matrix}$
where, as before, the reduced quantities with the subscript “y” and tilde above are defined using Ã_afrom (2.111) instead of A_a.
Even though K_Cin (2.100) can be any non-singular matrix, it is always possible to choose it in a way that the new state variables have particular physical significance. Similar to N_Lwhere the states can be transformed from currents i_xto fluxes λ_x, in the case of N_Cthe states may be transformed from voltages v_yto charges q_yby an appropriate choice of transformation. In particular, by defining the transformation as
K_C=C_y ⁻¹ (2.112)
v_y=C_y ⁻¹q_y (2.113)
The state equations (2.111) with capacitor charges being states becomes
$\begin{matrix} \frac{\partial q_{y}}{\partial t} = - G_{y} C_{y}^{- 1} q_{y} + A_{a} j_{br} & (2.114) \end{matrix}$

Multiple Networks

An algorithm for automatically deriving a system (or systems) of DAEs for time-domain simulation of many practically useful electrical networks will now be discussed. As noted previously, there are some conditions under which a global network with state variables, here denoted as
N=(G,P _ζ), where ζ=L, C (3.1)
cannot be structurally represented entirely only as N_Lor as N_C. Specifically, an attempt to represent the network N structurally entirely as N_Lfails if there is a shortage of inductive link branches in the set B_x, On the other hand, the network N cannot be viewed only as N_Cif the corresponding forest of spanning trees lacks capacitive branches in the set B_y. In these cases, the global network cannot be viewed entirely as a single type of a network, and therefore, a more general approach is required.
A method of handling networks with arbitrary topologies is presented here based on the separation of N into interconnected networks. That is, if it is not possible to represent an entire circuit as a network of a single kind, it is necessary to identify some cuts of N that can be grouped together to form several networks based on their topological properties. Separation of N into N_LN_C, and N_Ain a way such that it is possible to obtain consistent state equations for N_L, N_C, and a system of algebraic equations for N_Ais, therefore, a generalization of the ASMG approach for circuits with arbitrary topologies. In this sense, a robust state selection algorithm is a key to such a generalization.
The selection of branches that can carry state variables in N_Lor N_Cis determined through a topological partitioning of the global graph G^g=(N,B) into a forest of trees and remaining link branches as
G ^g=(N,B)={G _trees ,G _links} (3.2)
where
G _trees=(N,B _y) (3.3)
G _links=(N,B _x) (3.4)
In previous discussions devoted to networks with state variables, the actual partitioning of G^g=(N,B) was performed with two distinct objectives. First, with the help of the MinSTA and the weight function w_L(b), it was possible to reorder branches such the that set of link-branches B_xcould be chosen to represent state variables in N_L. Second, using a similar topological approach based on MaxSTA and the weight function w_C(b), it was shown that branches from the set of tree-branches B_ycan represent state variables in N_C. Thereafter, the conditions for the state equations to be complete were: for N_Lthe corresponding set B_xmust contain only inductive branches; and for N_Cthe corresponding set B_ymust contain only capacitive branches. These two branch sets could be obtained independently running MinST and MaxST algorithms with different weight functions. These topological conditions are stated in (2.82) and (2.87).
The technique of network identification and partitioning introduced here is based on the TCF of the node incidence matrix assembled for the global network. As it will be shown, the TCF makes topological information about N available in an “algebraic” sense for the further use in KCLs and KVLs. There are also some details about notation that are worth pointing out. The topological quantities referred to a particular branch order, such as in (2.17), will be distinguished by the hat sign, keeping in mind that it is always possible to transform them back to the original order using a corresponding permutation matrix. The superscripts “L”, “C”, and “A”, and the combinations thereof would be employed to relate variables and quantities to the inductive, capacitive, algebraic, and the overlapping networks, respectively.
Starting with N_L, in will be shown how an algebraic network can be identified and how the corresponding system of algebraic equations can be assembled. It will also be shown that even in the presence of N_A, a minimal and consistent state equation for N_Lcan still be obtained following the same dimensionality reduction procedure set forth in Wasynczuk and Sudhoff. Then, similar derivations will be repeated for N_C, in a somewhat simplified form, using its structural duality with respect to the inductive case. Finally, a way of obtaining a consistent system of DAEs relating all networks will be presented using the established framework.

Inductive and Algebraic Network Interconnection

Here, it is assumed that a network N=(G, P_L) is constructed using branches shown in FIGS. 4A and 4C. This time, it is assumed that the MinSTA with weight function w_L(b) is applied, and that in the end (2.82) does not hold. This implies that there is no way the set of links B_xcan be chosen to contain only inductive branches, and
$\begin{matrix} \sum_{j \in B_{x}} w_{L} (b_{j}) = h < size (B_{x}) & (3.5) \end{matrix}$
It follows from (2.17) and (3.5), there is at least one non-inductive branch in B_xaddition of which to G_trees=(N,B_y) would result in a cycle. This cycle would be entirely composed of non-inductive branches. Otherwise, the spanning tree corresponding to this link branch would not be minimal which, in turn, would contradict the result of the MinSTA. In general, addition of any branch from G_links=(N,B_x,) to G_trees=(N,B_y) results in a cycle. Therefore, the global graph G^g=(N, B) has as many non-inductive cycles as there are non-inductive branches in the branch set B_x. This property is based on the MinST and the weight function (2.13). Thus, based on (3.2), the total number of non-inductive cycles in G^g=(N,B) can be determined as
m=size(B _x)−h=s−h (3.6)
These cycles do not have state variables such as currents or fluxes. Therefore, none of the branches participating in such cycles can be placed in N_L.
There are many methods that can be used to identify cycle(s) of a particular kind in a graph. Recall that the branches can be reordered similar to (2.17) and (2.18) by a corresponding permutation matrix obtained from the MinSTA. This time, it is necessary to reorder branches in (2.17) and columns in (2.18) in an even more sophisticated way. Specifically, since it is known which branches in B_xare non-inductive, they can be identified and put first on the left side in A_links. Then, the branches in B_ythat correspond to cycles linked by non-inductive branches in B_xare identified, and the block A_treesis resorted such that these columns appear on the right. Similar to (2.18) this final branch ordering can be expressed in terms of the node incidence matrix and permutation matrix as
A_fT_L=[A_trees ^L,A_trees ^A,A_links ^A,A_links ^L] (3.7)
According to the system of notation employed, the superscript “L” denotes that the corresponding branches can be safely placed into inductive network N_L, and the “A” identifies all other non-inductive branches, as viewed from N_L, that belong to N_A. Taking the RREF of (3.7) and removing zero rows on the bottom yields a TCF similar to the one introduced in (2.21).
$\begin{matrix} A_{TCF} = [\begin{matrix} I_{η \times η} & 0 & 0 & {\hat{A}}_{η \times h}^{L} \\ 0 & I_{μ \times μ} & {\hat{A}}_{μ \times m}^{A} & {\hat{C}}_{μ \times h}^{A} \end{matrix}] & (3.8) \end{matrix}$
where
h—is the number of inductive link-branches in the set B_x, as defined in (3.5);
m—is the number of non-inductive link-branches in B_x, as defined in (3.6);
η—is the number of tree-branches in the set B_ythat are linked by the h inductive link-branches in B_x.
μ—is the number of tree-branches in B_ythat are linked by the m non-inductive link-branches from B_x.
Also, η+μ=r, and m+h=s, where s and r are defined in (2.8) and (2.9), respectively. The corresponding subsets of tree-branches in G_trees=(N,B_y) can be readily identified from the node incidence matrix A_fin its TCF such as (2.21).
The reduced node incidence matrix and the basic loop matrix for the global network N are found as usual
$\begin{matrix} {\hat{A}}_{a} = A_{TCF} & (3.9) \\ {\hat{B}}_{b} = \langle \begin{matrix} 0_{m \times η} & - {({\hat{A}}_{μ \times m}^{A})}^{T} & I_{m \times m} & 0_{m \times h} \\ - {({\hat{A}}_{η \times h}^{L})}^{T} & - {({\hat{C}}_{μ \times h}^{A})}^{T} & 0_{h \times m} & I_{h \times h} \end{matrix} \rangle & (3.10) \end{matrix}$
Even though the TCFs (2.21) and (3.8) are different only in terms of the relative order of some columns (branches), expression (3.8) has more advanced structural properties that can be utilized in writing KCL and KVL for the sections of the global network. Specifically, KCL for N_Lcan be written based on (3.9) as
$\begin{matrix} \begin{matrix} [I_{η \times η} 0_{η \times μ + m} {\hat{A}}_{η \times h}^{L}] {\hat{i}}_{br} = [I_{η \times η} {\hat{A}}_{r - μ \times h}^{L}] {\hat{i}}_{br}^{L} \\ = {\hat{A}}_{a}^{L} {\hat{i}}_{br}^{L} \\ = 0 \end{matrix} & (3.11) \end{matrix}$
Based on (3.10), KVL for N_Acan be derived as
$\begin{matrix} \begin{matrix} [0_{m \times η} - {({\hat{A}}_{μ \times m}^{A})}^{T} I_{m \times m} 0_{m \times h}] {\hat{v}}_{br} = [- {({\hat{A}}_{μ \times m}^{A})}^{T} I_{m \times m}] {\hat{v}}_{br}^{A} \\ = {\hat{B}}_{b}^{A} {\hat{v}}_{br}^{A} \\ = 0 \end{matrix} & (3.12) \end{matrix}$
It is interesting to observe that KCL (3.11) is a self-contained equation for N_L. Similarly, KVL (3.12) contains only quantities relevant to N_A. In this sense, these two equations are de-coupled.
The vector of state variables for N_L, which is the vector of independent currents, is then chosen as
$\begin{matrix} {\hat{i}}_{br}^{L} = [\begin{matrix} - {\hat{A}}_{η \times h}^{L} \\ I_{h \times h} \end{matrix}] i_{x} = {({\hat{B}}_{b}^{L})}^{T} i_{x} & (3.13) \end{matrix}$
In proceeding further, two more equations describing both networks are needed, specifically KVL for N_L, and KCL for N_A. First, from (3.10), the following KVL can be written
$\begin{matrix} ⌈ - {({\hat{A}}_{η \times h}^{L})}^{T} - {({\hat{C}}_{μ \times h}^{A})}^{T} 0_{h \times m} I_{h \times h} ⌉ {\hat{v}}_{br} = 0 & (3.14) \end{matrix}$
which may be rewritten using notation for the two networks as
$\begin{matrix} [- {({\hat{A}}_{η \times h}^{L})}^{T} I_{h \times h}] {\hat{v}}_{br}^{L} = [{({\hat{C}}_{μ \times h}^{A})}^{T} 0_{h \times m}] {\hat{v}}_{br}^{A} {\hat{B}}_{b}^{L} {\hat{v}}_{br}^{L} = {\hat{C}}_{b}^{LA} {\hat{v}}_{br}^{A} & (3.15) \end{matrix}$
Similarly, from (3.9), KCL can be written as
$\begin{matrix} [0_{μ \times η} I_{μ \times μ} {\hat{A}}_{μ \times m}^{A} {\hat{C}}_{μ \times h}^{A}] {\hat{i}}_{br} = 0 and & (3.16) \\ [I_{μ \times μ} {\hat{A}}_{μ \times m}^{A}] {\hat{i}}_{br}^{A} = - [0_{μ \times η} {\hat{C}}_{μ \times h}^{A}] {\hat{i}}_{br}^{L} {\hat{A}}_{a}^{A} {\hat{i}}_{br}^{A} = {\hat{C}}_{a}^{LA} {\hat{i}}_{br}^{L} = - {\hat{C}}^{A} i_{x} & (3.17) \end{matrix}$
It is worth noting that (3.15) and (3.17) express the connection between the two networks. To be more specific, KVL (3.15) for N_Lnow has an additional voltage source term that comes from the interconnection of N_Land N_A. Similarly, KCL (3.17) written for N has current source term which comes from N_L.
The dimensionality reduction procedure for obtaining state equations for N_Lremains the same. That is, the voltage equation written in matrix form for the N is multiplied from the left by the corresponding KVL loop matrix, and the vector of branch currents is substituted with an expression of the form (2.83). Thus, using B_b ^Lfrom the KVL (3.15) and the vector of state currents chosen as in (3.13), the state equation for N_Lbecomes
$\begin{matrix} \frac{\partial i_{x}}{\partial t} = - L_{x}^{- 1} (R_{x} + \frac{\partial L_{x}}{\partial t}) i_{x} - L_{x}^{- 1} B_{b}^{L} e_{br}^{L} - L_{x}^{- 1} C_{b}^{LA} v_{br}^{A} & (3.18) \end{matrix}$
where the reduced quantities with the subscript “x” can be obtained as
$\begin{matrix} \begin{matrix} R_{x} = B_{b}^{L} {R_{br}^{L} (B_{b}^{L})}^{T} \\ = [- {({\hat{A}}_{η \times h}^{L})}^{T} \begin{matrix} 0_{h \times μ + m} & I_{h \times h} \end{matrix}] {{\hat{R}}_{br} [- {({\hat{A}}_{η \times h}^{L})}^{T} \begin{matrix} 0_{h \times μ + m} & I_{h \times h} \end{matrix}]}^{T} \end{matrix} & (3.19) \\ \begin{matrix} L_{x} = B_{b}^{L} {L_{br}^{L} (B_{b}^{L})}^{T} \\ = [- {({\hat{A}}_{η \times h}^{L})}^{T} \begin{matrix} 0_{h \times μ + m} & I_{h \times h} \end{matrix}] {{\hat{L}}_{br} [- {({\hat{A}}_{η \times h}^{L})}^{T^{'}} \begin{matrix} 0_{h \times μ + m} & I_{h \times h} \end{matrix}]}^{T} \end{matrix} & (3.20) \\ R_{br} = T_{L}^{T} R_{br} T_{L} & (3.21) \\ L_{br} = T_{L}^{T} L_{br} T_{L} & (3.22) \end{matrix}$
It can be verified that L_xis non-singular, and that (3.18) is in fact a minimal state equation for N_Lwith the total of h state variables. Also, note the additional source term on the right side of (3.18), which distinguishes (3.18) from (2.84).
It appears that in order to solve (3.18), the branch voltages for N_Ashould be known. These voltages are functions of internal topology, parameters, and external sources. In general, it is necessary to compute both currents and voltages for all branches in N_A. Here, KVL (3.12) and KCL (3.17) are utilized. In particular, suppose the currents of N_Afirst, and then the branch voltages. However, in (3.17) there are fewer equations than branches in N_A. Specifically, (3.17) provides μ equations with m+μ unknowns. That is, there should be m more equations, precisely one for each cycle in N_A. On the other hand, N_Ahas its own voltage equation, which may be expressed in terms of the branch order given by the columns of (3.7) as
{circumflex over (v)} _br ^A ={circumflex over (R)} _br ^A î _br ^A +ê _br ^A (3.23)
Substituting (3.23) into (3.11), the following result can be obtained
{circumflex over (B)} _b ^A {circumflex over (R)} _br ^A î _br ^A=−{circumflex over (B)}_b ^A ê _br ^A (3.24)
where ê_br ^Ais the vector of voltage sources corresponding to branches in N_Aand the branch order given by (3.7).
Combining (3.17) and (3.24), a complete system of m+μ equations for the branch currents of N_Acan be assembled in the following form
$\begin{matrix} [\begin{matrix} {\hat{A}}_{a}^{A} \\ {\hat{B}}_{b}^{A} {\hat{R}}_{br}^{A} \end{matrix}] {\hat{i}}_{br}^{A} = - [\begin{matrix} {\hat{C}}^{A} & 0 \\ 0 & {\hat{B}}_{b}^{A} \end{matrix}] [\begin{matrix} i_{x} \\ {\hat{e}}_{br}^{A} \end{matrix}] & (3.25) \end{matrix}$
If the branch resistances do not depend on currents, the system of equations (3.25) has the form Ax=b. In addition, if N_Ahas time-invariant resistive branches, A is a constant nonsingular matrix. Therefore, for the time-invariant case, (3.25) can be solved by inverting the corresponding matrix A once for each new topology. Thereafter, the corresponding vector of branch voltages {circumflex over (v)}_br ^Ais found using (3.23). Finally, (3.18), (3.23), and (3.25) form a system of DAEs describing the two networks N_Land N_A.

Capacitive and Algebraic Network Interconnection

In this section, it is assumed that a network N=(G, P_C) is constructed from branch models shown in FIG. 4. It is assumed that the MaxSTA with the weight function w_C(b) is applied and that condition (2.87) does not hold. This result implies that the set of tree-branches B_ycannot be chosen to contain only capacitive branches, and
$\begin{matrix} \sum_{i \in B_{y}} w_{C} (b_{i}) = τ < size (B_{y}) & (3.26) \end{matrix}$
Then, based on (2.17) and (3.26), the total number of non-capacitive tree-branches in the forest of maximum spanning trees can be determined as
ζ=size(B _y)−τ=r−τ (3.27)
When applying the MaxSTA, a permutation matrix that sorts columns of the full node incidence matrix A_fsimilar to (2.18) is also assembled. Using (2.18) and the maximum spanning tree property, the branches in A_linksand A_treesare reordered in the following way.
First, all τ columns of A_links, corresponding to capacitive branches are placed on the left of A_links. Then, from the set B_x, the branches that are links to the capacitive trees in G_trees=(N,B_x) are identified and the corresponding columns put on the right of A_trees. The final branch ordering with the corresponding permutation matrix yields a result similar to (3.7). That is
A_fT_C[A_trees ^C,A_trees ^A,A_links ^A,A_links ^C] (3.28)
Another TCF can be obtained by taking the RREF of (3.28) and removing zero rows from the bottom. This TCF has the following structure
$\begin{matrix} A_{TCF} = [\begin{matrix} I_{τ \times τ} & 0 & {\hat{D}}_{τ \times t}^{CA} & {\hat{A}}_{τ \times z}^{C} \\ 0 & I_{ζ \times ζ} & {\hat{A}}_{ζ \times t}^{A} & 0 \end{matrix}] & (3.29) \end{matrix}$
where
t—is the number of non-capacitive link-branches in the set B_x.
z—is the number of capacitive link-branches in the set B_x.
τ—is the number of capacitive tree-branches in the set B_yas defined in (3.26).
ζ—is the number of non-capacitive tree-branches in the set B_yas defined in (3.27).
Also, ζ+τ=r and t+z=s. If any of the dimensions are equal to zero, the corresponding block in (3.29) would simply disappear.
Based on the TCF (3.29), the reduced node incidence matrix and the basic loop matrix, both referred to the branch ordering defined in (3.28), are found as
$\begin{matrix} A_{a} = A_{TCF} & (3.30) \\ {\hat{B}}_{b} = \langle \begin{matrix} - {({\hat{D}}_{τ \times t}^{CA})}^{T} & - {({\hat{A}}_{ζ \times t}^{A})}^{T} & I_{t \times τ} & 0 \\ \leftarrow {({\hat{A}}_{τ \times z}^{C})}^{T} & 0 & 0 & I_{z \times z} \end{matrix} \rangle & (3.31) \end{matrix}$
Similar to (3.11)-(3.17), KCLs and KVLs can be written for the two networks. That is, writing KCL using (3.30),
$\begin{matrix} [\begin{matrix} I_{τ \times τ} & 0 & \begin{matrix} {\hat{D}}_{τ \times t}^{CA} & {\hat{A}}_{τ \times z}^{C} \end{matrix} \end{matrix}] {\hat{i}}_{br} = 0 & (3.32) \\ [\begin{matrix} I_{τ \times τ} & {\hat{A}}_{τ \times z}^{C} \end{matrix}] {\hat{i}}_{br}^{C} = [\begin{matrix} 0 & - {\hat{D}}_{τ \times t}^{CA} \end{matrix}] {\hat{i}}_{br}^{A} {\hat{A}}_{a}^{C} {\hat{i}}_{br}^{C} = {\hat{D}}_{a}^{CA} {\hat{i}}_{br}^{A} and & (3.33) \\ [I_{ζ \times ζ} {\hat{A}}_{ζ \times t}^{A}] {\hat{i}}_{br}^{A} = {\hat{A}}_{a}^{A} {\hat{i}}_{br}^{A} = 0 & (3.34) \end{matrix}$
Writing KVL based on (3.31),
$\begin{matrix} ⌈ - {({\hat{D}}_{τ \times t}^{A})}^{T} - {({\hat{A}}_{ζ \times t}^{A})}^{T} \begin{matrix} I_{t \times t} & 0_{t \times z} \end{matrix} ⌉ {\hat{v}}_{br} = 0 & (3.35) \\ [- {({\hat{A}}_{ζ \times t}^{A})}^{T} I_{t \times t}] {\hat{v}}_{br}^{A} = [{({\hat{D}}_{τ \times t}^{A})}^{T} 0_{z \times z}] {\hat{v}}_{br}^{C} {\hat{B}}_{b}^{A} {\hat{v}}_{br}^{A} = {\hat{D}}_{b}^{CA} {\hat{v}}_{br}^{C} and & (3.36) \\ ⌈ - {({\hat{A}}_{τ \times z}^{C})}^{T} I_{z \times z} ⌉ {\hat{v}}_{br}^{C} = {\hat{B}}_{b}^{C} {\hat{v}}_{br}^{C} = 0 & (3.37) \end{matrix}$
Selecting the vector of state variables for N_Cto be the vector of independent voltages v_y, similar to (2.88), then
$\begin{matrix} {\hat{v}}_{br}^{C} = ⌈ \begin{matrix} I_{τ \times τ} \\ {({\hat{A}}_{τ \times z}^{C})}^{T} \end{matrix} ⌉ v_{y} = {({\hat{A}}_{a}^{C})}^{T} v_{y} & (3.38) \end{matrix}$
Expression (3.28) can also be written as
{circumflex over (B)} _b ^A {circumflex over (v)} _br ^A=({circumflex over (D)} ^A)^T v _y (3.39)
Assuming that the vector of branch currents for network N_Acan be expressed in terms of the branch order (3.28) as
î _br ^A =Ĝ _br ^A {circumflex over (v)} _br ^A −ĵ _br ^A (3.40)
the system of linear equations for {circumflex over (v)}_brbased on (3.34), (3.39), and (3.40) can be written as
$\begin{matrix} [\begin{matrix} {\hat{B}}_{b}^{A} \\ {\hat{A}}_{a}^{A} {\hat{G}}_{br}^{A} \end{matrix}] {\hat{v}}_{br}^{A} = ⌈ \begin{matrix} {({\hat{D}}^{A})}^{T} & 0 \\ 0 & {\hat{A}}_{a}^{A} \end{matrix} ⌉ [\begin{matrix} v_{y} \\ {\hat{j}}_{br}^{A} \end{matrix}] & (3.41) \end{matrix}$
If the branch conductances do not depend on voltages, the system of equations (3.41) has the form Ax=b, with ζ+t equations and the same number of unknowns, and A being a full-rank square matrix.
Using the dimensionality reduction procedure in Wasynczuk and Sudhoff, the state equation for the N_Cin the presence of the algebraic network becomes
$\begin{matrix} \frac{\partial v_{y}}{\partial t} = - C_{y}^{- 1} (G_{y} + \frac{\partial C_{y}}{\partial t}) v_{y} + C_{y}^{- 1} A_{a} j_{br} - C_{y}^{- 1} D_{a}^{CA} i_{br}^{A} & (3.42) \end{matrix}$
where the reduced quantities are defined as
$\begin{matrix} \begin{matrix} C_{y} = A_{a}^{C} {C_{br} (A_{a}^{C})}^{T} \\ = [I_{τ \times τ} 0_{τ \times ζ + t} {\hat{A}}_{τ \times z}^{C}] {{\hat{C}}_{br} [I_{τ \times τ} 0_{τ \times ζ + t} {\hat{A}}_{τ \times z}^{C}]}^{T} \end{matrix} & (3.43) \\ \begin{matrix} G_{y} = A_{a}^{C} {G_{br} (A_{a}^{C})}^{T} \\ = [I_{τ \times τ} 0_{τ \times ζ + t} {\hat{A}}_{τ \times z}^{C}] {{\hat{G}}_{br} [I_{τ \times τ} 0_{τ \times ζ + t} {\hat{A}}_{τ \times z}^{C}]}^{T} \end{matrix} & (3.44) \end{matrix}$
and
Ĉ_br=T_L ^TC_brT_L (3.45)
Ĝ_br=T_L ^TG_brT_L (3.46)

Inductive, Capacitive, and Algebraic Network Interconnection

Heretofore, both types of network, namely N_Land N_C, have been considered with a shortage of branches that can carry state variables. In both cases, the corresponding minimal state equations (3.18) and (3.42) were completed by adding extra source terms that came about due to the algebraic part of the corresponding network. Also, both types of the final systems of DAEs are structurally similar due to their dual circuit. These structural properties can be utilized even further in obtaining a system of DAEs for the global network N=(G, P) that includes N_L, N_C, and N_A. In doing so, there are two approaches.
In the first approach, N_Land N_Aare considered as was done in the beginning of this chapter, whereupon N_Cis incorporated into the existing structure. In doing so, capacitor voltages of N_Ccan be mapped into e_br ^Land v_br ^Ain (3.18). A second approach would consist of adding N_Lto the structure developed for the N_Cand N_A, and mapping inductor currents into j_br ^Cand i_br ^Ain (3.42). Since both methods yield equivalent results, choosing either one is a matter of pure convenience. Following the order in which the material was presented earlier, preference is given to the first approach.
Suppose that the MinSTA is applied resulting in (3.5) from which the branch order (3.7) and the TCF (3.8) are obtained. In the most general case, N_Ccould have its branches anywhere in G_trees=(N,B_y) and among any non-inductive link-branches in G_links=(N,B_x). Relative to (3.7), these branches could be columns corresponding to blocks A_trees ^L, A_trees ^A, and A_links ^A. The capacitive branches corresponding to any of the columns in A_trees ^Lrepresent the part of N_Cthat overlaps with N_L. As a result of such overlapping, each of such capacitive branches is going to have its own independent state variable within N_C. In fact, the branch voltages corresponding to such capacitive branches can be viewed as independent voltage sources e_br ^Lpresent in (3.18). The remaining capacitive branches are viewed as a part of an algebraic network for N_L, and therefore, are going to be represented in the columns of blocks A_trees ^Aand A_links ^Aas given in (3.7).
Now, the challenge is to reorder columns (branches) in (3.7) taking into consideration the capacitive network. The procedure of reordering columns is very similar to the two previous cases. In particular, all columns corresponding to the non-capacitive branches in A_trees ^Lare identified and placed on the left side of this block. Then, the trees of capacitive branches that form a capacitive network need to be separated. In order to achieve this, the MaxSTA with weight function w_C(b) is applied to the branches in A_trees ^Aand A_links ^A. Then, with the result (3.26), the columns of this block are again sorted similar to (3.28). Thereafter, the final branch ordering with corresponding permutation matrix may be expressed as
A_fT_LCA=[A_trees ^L,A_trees ^LC,A_trees ^C,A_trees ^A,A_links ^A,A_links ^CA,A_links ^L] (3.47)
The permutation matrix T_LCAin (3.47) sorts branches of the global network N in groups with very specific topological properties corresponding to different networks. Again, taking the RREF of (3.47) and removing the zero rows, a TCF with the following structure is produced
$\begin{matrix} A_{TCF} = [\begin{matrix} I_{η \times η}^{L} & 0 & 0 & 0 & 0 & 0 & {\hat{A}}_{η \times h}^{L} \\ 0 & I_{μ \times μ}^{LC} & 0 & 0 & 0 & 0 & {\hat{A}}_{μ \times h}^{LC} \\ 0 & 0 & I_{τ \times τ}^{CA} & 0 & {\hat{D}}_{τ \times t}^{CA} & {\hat{A}}_{τ \times z}^{C} & {\hat{C}}_{τ \times h}^{LC} \\ 0 & 0 & 0 & I_{ζ \times ζ}^{A} & {\hat{A}}_{ζ \times t}^{A} & 0 & {\hat{C}}_{ζ \times h}^{LA} \end{matrix}] & (3.48) \end{matrix}$
where
η—is the number of non-capacitive tree-branches in the set B_ythat are also placed in N_L.
μ—is the number of capacitive tree-branches in B_ythat are placed in N_Land N_C.
τ—is the number of capacitive tree-branches in B_ythat are placed in N_C.
ζ—is the number of non-capacitive tree-branches in B_ythat are placed in N_A.
t—is the number of non-capacitive link branches in B_x, that are placed in N_A.
z—is the number of capacitive link-branches in B_x, that are placed in N_C.
h—is the number of inductive link-branches in the set B_xthat are placed in N_L.
The reduced node incidence matrix and the basic loop matrix for the global network referred to the branch order (3.47) are then found as
$\begin{matrix} A_{a} = A_{TCF} & (3.49) \\ {\hat{B}}_{b} = [\begin{matrix} 0 & 0 & - {({\hat{D}}_{τ \times t}^{CA})}^{T} & - {({\hat{A}}_{ζ \times t}^{A})}^{T} & I_{t \times t} & 0 & 0 \\ 0 & 0 & - {({\hat{A}}_{τ \times z}^{C})}^{T} & 0 & 0 & I_{z \times z} & 0 \\ - {({\hat{A}}_{η \times h}^{L})}^{T} & - {({\hat{A}}_{μ \times h}^{LC})}^{T} & - {({\hat{C}}_{τ \times h}^{LC})}^{T} & - {({\hat{C}}_{ζ \times h}^{LA})}^{T} & 0 & 0 & I_{h \times h} \end{matrix}] & (3.50) \end{matrix}$
From this point, it is possible to proceed as usual. KCL for N_Lusing (3.49) can be written as
$\begin{matrix} [\begin{matrix} I_{η \times η}^{L} & 0 & 0_{η \times τ + ζ + t + z} & {\hat{A}}_{η \times h}^{L} \\ 0 & I_{μ \times μ}^{LC} & 0_{μ \times τ + ζ + t + z} & {\hat{A}}_{μ \times h}^{LC} \end{matrix}] {\hat{i}}_{br} = 0 & (3.51) \\ [\begin{matrix} I_{η \times η}^{L} & 0 & {\hat{A}}_{η \times h}^{L} \\ 0 & I_{μ \times μ}^{LC} & {\hat{A}}_{μ \times h}^{LC} \end{matrix}] {\hat{i}}_{br}^{L} = {\hat{A}}_{a}^{L} {\hat{i}}_{br}^{L} = 0 & (3.52) \end{matrix}$
For N_Land N_A,
$\begin{matrix} [0_{ζ \times η + μ + τ} I_{ζ \times ζ}^{A} {\hat{A}}_{ζ \times t}^{A} 0_{ζ \times z} {\hat{C}}_{ζ \times h}^{LA}] {\hat{i}}_{br} = 0 & (3.53) \\ [I_{ζ \times ζ}^{A} {\hat{A}}_{ζ}^{A} \times t] {\hat{i}}_{br}^{A} = [0_{ζ \times η + μ + τ} - {\hat{C}}_{ζ \times h}^{LA}] {\hat{i}}_{br}^{L} {\hat{A}}_{a}^{A} {\hat{i}}_{br}^{A} = - {\hat{C}}_{a}^{LA} {\hat{i}}_{br}^{L} & (3, 54) \end{matrix}$
Finally, KCL relating N_Cto the two remaining networks is
$\begin{matrix} [\begin{matrix} 0_{μ \times η} & I_{μ \times μ}^{LC} & 0 & 0 & 0 & 0 & {\hat{A}}_{μ \times h}^{LC} \\ 0_{τ \times η} & 0_{τ \times μ} & I_{τ \times τ}^{CA} & 0_{τ \times ζ} & {\hat{D}}_{τ \times t}^{CA} & {\hat{A}}_{τ \times z}^{C} & {\hat{C}}_{τ \times h}^{LC} \end{matrix}] {\hat{i}}_{br} = 0 & (3.55) \\ [\begin{matrix} I_{μ \times μ}^{LC} & 0 & 0 \\ 0_{τ \times μ} & I_{τ \times τ}^{CA} & {\hat{A}}_{τ \times z}^{C} \end{matrix}] {\hat{i}}_{br}^{C} = [\begin{matrix} 0_{μ \times η + μ} & - {\hat{A}}_{μ \times h}^{LC} \\ 0_{τ \times η + μ} & - {\hat{C}}_{τ \times h}^{LC} \end{matrix}] {\hat{i}}_{br}^{L} + [\begin{matrix} 0_{μ \times ζ} & 0_{μ \times t} \\ 0_{τ \times ζ} & - D_{τ \times t} \end{matrix}] {\hat{i}}_{br}^{A} {\hat{A}}_{a}^{C} {\hat{i}}_{br}^{C} = - {\hat{C}}_{a}^{LC} {\hat{i}}_{br}^{L} - {\hat{D}}_{a}^{CA} {\hat{i}}_{br}^{A} & (3.56) \end{matrix}$
Based on (3.50), KVL can be written in the following way. First, for the capacitive network
┌0_z×η+μ−(Â _τ×z ^C)^T0_z×τ+t I _z×z0_z×h ┐{circumflex over (v)} _br=0 (3.57)
┌−(Â _τ×z ^C)^T I _z×z ┐{circumflex over (v)} _br ^C ={circumflex over (B)} _b ^C {circumflex over (v)} _br ^C=0 (3.58)

Then, for N_Cand N_A

$\begin{matrix} ⌈ 0_{τ \times η + μ} - {({\hat{D}}_{τ \times t}^{CA})}^{T} - {({\hat{A}}_{ζ \times t}^{A})}^{T} I_{t \times t} 0_{t \times z + h} ⌉ {\hat{v}}_{br} = 0 & (3.59) \\ [- {({\hat{A}}_{ζ \times t}^{A})}^{T} I_{t \times t}] {\hat{v}}_{br}^{A} = [0_{t \times μ} {({\hat{D}}_{τ \times t}^{CA})}^{T} 0_{t \times z}] {\hat{v}}_{br}^{C} {\hat{B}}_{b}^{A} {\hat{v}}_{br}^{A} = - {\hat{D}}_{b}^{CA} {\hat{v}}_{br}^{C} & (3.60) \end{matrix}$
Finally, KVL relating all three networks can be written as
$\begin{matrix} ⌈ - {({\hat{A}}_{η \times h}^{L})}^{T} - {({\hat{A}}_{μ \times h}^{LC})}^{T} - {({\hat{C}}_{τ \times h}^{LC})}^{T} - {({\hat{C}}_{ζ \times h}^{LA})}^{T} 0_{h \times t + z} I_{h \times h} ⌉ {\hat{v}}_{br} = 0 & (3.61) \\ \begin{matrix} [- {({\hat{A}}_{η \times h}^{L})}^{T} 0_{h \times μ} I_{h \times h}] {\hat{v}}_{br}^{L} = [{({\hat{A}}_{μ \times h}^{LC})}^{T} {({\hat{C}}_{τ \times h}^{LC})}^{T} 0_{h \times z}] {\hat{v}}_{br}^{C} + ⌈ {({\hat{C}}_{ζ \times h}^{LA})}^{T} 0_{2} \dots ⌉ {\hat{v}}_{br}^{A} \\ = {\hat{B}}_{b}^{L} {\hat{v}}_{br}^{L} \\ = - {\hat{C}}_{b}^{LC} {\hat{v}}_{br}^{C} - {\hat{C}}_{b}^{LA} {\hat{v}}_{br}^{A} \end{matrix} & (3.62) \end{matrix}$
As before, the vector of state variables for N_Lis selected as a vector of independent currents such that
$\begin{matrix} {\hat{i}}_{br}^{L} = [\begin{matrix} - {\hat{A}}_{η \times h}^{L} \\ - {\hat{A}}_{μ \times h}^{LC} \\ I_{h \times h} \end{matrix}] i_{x} = {({\hat{B}}_{b}^{LC})}^{T} i_{x} & (3.63) \end{matrix}$
Note that in (3.63), it is possible to use {circumflex over (B)}_b ^Lfrom (3.62) in order to avoid computing currents of capacitive branches.
Similar 016410-112 to (3.38), the vector of states for N_Cis chosen to be a vector of independent capacitor voltages such that
$\begin{matrix} {\hat{v}}_{br}^{C} = ⌈ \begin{matrix} I_{μ \times μ}^{LC} & 0 \\ 0 & I_{τ \times τ}^{CA} \\ 0 & {({\hat{A}}_{τ \times z}^{C})}^{T} \end{matrix} ⌉ v_{y} = {({\hat{A}}_{a}^{C})}^{T} v_{y} & (3.64) \end{matrix}$
Based on definition (3.63), expression (3.54) becomes
$\begin{matrix} {\ddot{A}}_{a}^{A} {\hat{i}}_{br}^{A} = - {\ddot{C}}_{a}^{LA} {\hat{i}}_{br}^{L} = - C^{LA} i_{x} & (3.65) \end{matrix}$
Also, based on (3.63) and (3.56), it follows that
$\begin{matrix} [\begin{matrix} 0_{μ \times η + μ} & - {\hat{A}}_{μ \times h}^{LC} \\ 0_{τ \times η + μ} & - {\hat{C}}_{τ \times h}^{LC} \end{matrix}] {\hat{i}}_{br}^{L} = [\begin{matrix} - {\hat{A}}_{μ \times h}^{LC} \\ - {\hat{C}}_{τ \times h}^{LC} \end{matrix}] i_{x} = - D^{LC} i_{x} & (3.66) \end{matrix}$
which then can be used to rewrite the compact form of (3.56) as
$\begin{matrix} {\hat{A}}_{a}^{C} {\hat{i}}_{br}^{C} = - {\hat{C}}_{a}^{LC} {\hat{i}}^{L_{br}} - {\hat{D}}_{a}^{CA} {\hat{i}}_{br}^{A} = - D^{LC} i_{x} - {\hat{D}}_{a}^{CA} {\hat{i}}_{br}^{A} & (3.67) \end{matrix}$
Similarly, based on (3.64), KVL (3.60) and (3.62) can be rewritten as
$\begin{matrix} {\hat{B}}_{b}^{A} {\hat{v}}_{br}^{A} = - {\hat{D}}_{b}^{CA} {\hat{v}}_{br}^{C} = - D^{CA} v_{y} & (3.68) \\ ⌈ {({\hat{A}}_{μ \times h}^{LC})}^{T} {({\hat{C}}_{τ \times h}^{LC})}^{T} 0_{h \times z} ⌉ {\hat{v}}_{br}^{C} = ⌈ {({\hat{A}}_{μ \times h}^{LC})}^{T} {({\hat{C}}_{τ \times h}^{LC})}^{T} ⌉ v_{y} = - C^{LC} v_{y} & (3.69) \\ {\hat{B}}_{b}^{L} {\hat{v}}_{br}^{L} = - {\hat{C}}_{b}^{LC} {\hat{v}}_{br}^{C} - {\hat{C}}_{b}^{LA} {\hat{v}}_{br}^{A} = - C^{LC} v_{y} - {\hat{C}}_{b}^{LA} {\hat{v}}_{br}^{A} & (3.70) \end{matrix}$
Note that KCL (3.67) written for N_Cand KVL (3.70) written for N_Lare coupled through the corresponding interconnection matrices that are related to each other as
C ^LC=−(D ^LC)^T (3.71)
Using (3.70), (3.67), and the dimensionality reduction procedure of Wasynczuk and Sudhoff, the state equations for inductive and capacitive networks become
$\begin{matrix} \frac{\partial i_{x}}{\partial t} = - L_{x}^{- 1} (R_{x} + \frac{\partial L_{x}}{\partial t}) i_{x} - L_{x}^{- 1} C^{LC} v_{y} - L_{x}^{- 1} C_{b}^{LA} v_{br}^{A} - L_{x}^{- 1} B_{b}^{L} e_{br}^{L} & (3.72) \\ \frac{\partial v_{y}}{\partial t} = - C_{y}^{- 1} (G_{y} + \frac{\partial C_{y}}{\partial t}) v_{y} - C_{y}^{- 1} D^{LC} i_{x} - C_{y}^{- 1} D_{a}^{CA} i_{br}^{A} + C_{y}^{- 1} A_{a}^{C} j_{br}^{C} & (3.73) \end{matrix}$
where all reduced quantities with the subscripts “x” and “y” can be computed very much like the ones in (3.19), (3.20), (3.43), and (3.44) using the corresponding KVL matrix B_b ^Land KCL matrix A_a ^C. As before, it can be verified that (3.72) and (3.73) are the minimal state equations for the N_Land N_C, respectively. The reduced quantities L_xand C_yare indeed non-singular matrices for any practical network.
As before, in order to solve (3.72) and (3.73) for the state derivatives, the branch currents and voltages of the algebraic network must be computed first. Here, it is reasonable to assume that N_Acontains non-empty resistive branches that are modeled as depicted in FIGS. 4( a) and (b), with voltage or current sources, respectively.
Since a single branch cannot contain both types of sources and assuming that each branch of N_Ahas a nonzero resistor, the voltage equation for N_Acan be written as
{circumflex over (v)} _br ^A ={circumflex over (R)} _br ^A(î _br ^A +ĵ _br ^A)+ê _br ^A (3.74)
Then, in order to obtain a system of equations for N_A, KCL (3.65) and KVL (3.68) should be utilized. Thus, substituting (3.74) into (3.68) and combining the result with (3.65), the following system of equations is obtained
$\begin{matrix} [\begin{matrix} {\hat{A}}_{a}^{A} \\ {\hat{B}}_{b}^{A} {\hat{R}}_{br}^{A} \end{matrix}] {\hat{i}}_{br}^{A} = ⌈ \begin{matrix} - C^{LA} i_{x} \\ - D^{CA} v_{y} - {\hat{B}}_{b}^{A} ({\hat{R}}_{br}^{A} {\hat{j}}_{br}^{A} + {\hat{e}}_{br}^{A}) \end{matrix} ⌉ & (3.75) \end{matrix}$
which has ζ+t equations and the same number of unknowns. If {circumflex over (R)}_br ^Adoes not depend on branch currents, the system (3.75) would have the form Ax=b with A being a full rank square matrix. Thereafter, (3.75) can be solved for the vector of currents î_br ^A, which is then used in (3.74) in order to compute the voltages {circumflex over (v)}_br ^AFinally, (3.72)-(3.75) represent a consistent system of DAEs for the three elementary networks N_L, N_C, and N_A.

Elementary Network Identification Procedure

For large electrical networks, storing the associated multi-graphs in matrix format leads to sparse matrices of accordingly large dimensions. Even though the actual topological matrices used for elementary networks as KCLs and KVLs are smaller in size, these matrices are sparse and their dimensions are proportional to the size of the elementary networks. However, in order to gain computational efficiency, it is desirable to employ techniques that would avoid multiplication of large sparse matrices. A direct way to reduce the unnecessary operations is to use the techniques developed for sparse algebra in all computations involving topological matrices. Essentially this would lead to component-wise operations involving all of the network branch parameters. Alternatively, it is possible to use a more compact technique for graph representation such as collections of adjacency lists instead of adjacency matrices. Such lists can be conveniently implemented on a computer, which in turn would enable the usage, of more efficient Min/MaxST algorithms. Representation of graphs in terms of compact data structures also suggests reformulation of the network identification procedure using ordered sets instead of topological matrices.

Input Data

Given a general electrical network N, the associated graph is defined in terms of node and branch sets as G=(N,B). The graph G must be closed (G must form a valid circuit), which implies the following
(1) each branch in B is a member of some cycle composed of more than one branch (self-loops are not allowed)
(2) The entire branch set B can be partitioned into two subsets B_yand B_xsuch that they have no branches in common, set B_yhas a minimal size and spans the entire node set N. Moreover, G_trees=(N,B_y) is a spanning tree. If G is a multi-graph, then G_Tis a spanning forest (forest of spanning trees). Also, based on the parameter set P, the weight functions ω_L(2.81) and ω_C(2.86) must be defined over the entire branch set B. Then, depending on the order in which the elementary networks are to be found, a network identification procedure can be formulated in four major steps. Two procedures will now be discussed for the different order of network identification.

L-C-A Procedure

Step 1: Call the Minimum Spanning Tree Algorithm

MinSTA(G,w _L)=
G _trees=( N,B _y) (3.80)
Given B_yfrom (3.80), the set of link-branches is defined as
B _x =B−B _y (3.81)
Branches in B_xfor which ω_L(b_j)=1 (inductive link-branches) are identified defining the set B_x ^L. Then, the set of link-branches is re-ordered as follows
B_x={{tilde over (B)}_x ^CA,B_x ^L} (3.82)
where {circumflex over (B)}_x ^CAis a temporary set of non-inductive link-branches.
Step 2: B_x ^Lis removed from the global branch set B and the remaining closed graph is defined {tilde over (G)}=(Ñ, {tilde over (B)}^CA). Then, it is necessary to find all the tree-branches in B_ythat are linked by {tilde over (B)}_x ^CA. Let this temporary branch set be denoted as {tilde over (B)}_y ^CA. Thereafter, the branch set for G is found as
{tilde over (B)}^CA={tilde over (B)}_y ^CA,{circumflex over (B)}_x ^CA (3.83)
Based on the MinST property, it can be proved that none of the branches in {tilde over (B)}^CAare inductive. The branch set B_ycan be partitioned into the following sets
B_y={B_y ^LC,{tilde over (B)}_y ^CA} (3.84)
Since some of the branches in B_y ^LCmay be capacitive, this branch set may be partitioned further as
B_y ^LC={B_y ^L,B_y ^C ¹} (3.85)
Also note that since the branch set {tilde over (B)}^CA(3.83) forms a closed graph, the remaining branches {{tilde over (B)}_y ^LC,B_x ^L} do not, unless {tilde over (B)}^CAand {{tilde over (B)}_y ^LC, B_x ^L} are sets of galvanically disjointed branches.

Step 3: The Maximum Spanning Tree Algorithm is Applied

MaxSTA({tilde over (G)},w _C)=
{tilde over (G)} _trees=(Ñ,B _y ^CA) (3.86)
Since B_y ^CAhas an optimal property, it is no longer a temporary set. Also, due to such possible improvement, this set may be different from the one used in (3.83). Proceeding further, all branches in B_y ^CAfor which ω_C(b_j)=1 (capacitive branches) are identified and placed into a separate set here denoted as B_y ^C ²Then B_y ^CA; can be partitioned as
B_y ^CA={B_y ^C ²,B_y ^A} (3.87)
The set of link-branches is retrieved as
B _x ^CA ={tilde over (B)} ^CA −B _y ^CA (3.88)
Step 4: B_y ^Ais removed from the global branch set B_y ^CAand the remaining closed graph is denoted Ĝ=({circumflex over (N)},{circumflex over (B)}^C). Next, all link-branches in B_xthat correspond to tree-branches in B_y ^C ²are found. Let this branch set be denoted as B_x ^C. This set actually may or may not have capacitive branches, but all of its branches are links to capacitive trees in B_y ^C ². The remaining link-branches can be found as
B _x ^A =B _x ^CA −B _x ^C (3.89)
Then, the following set partitioning can be formed
B_y ^CA={B_y ^C ²,B_y ^A} (3.90)
B_x ^CA={B_x ^A,B_x ^C} (3.91)
Based on the property of MaxSTA, it is possible prove that B_x ^Ahas no capacitive branches.
Finally, the global branch set B can be reordered in terms of the sets identified so far as
B={B_y ^L, B_y ^C ¹, B_y ^C ², B_y ^A,B_x ^A, B_x ^C, B_x ^L} (3.92)
Thereafter, the elementary networks are formed based on the branch sets as follows
{B_y ^L, B_x ^L}∈N_L (3.93)
{B_y ^C ¹,B_y ^C ²,B_x ^C}∈N_C (3.94)
{B_y ^A,B_x ^A}∈N_A (3.95)
The branch order established in (3.92) corresponds to the TCF (3.48) with the block matrices having dimensions of their respective sets in (3.93)-(3.95).

C-L-A Procedure

Step 1: The Maximum Spanning Tree Algorithm is Applied

MaxSTA(G,w _C)=
G _trees=(N,B _y) (3.96)
The set of link-branches becomes
B _x =B−B _y (3.97)
All capacitive branches in B_yare identified. This set is then re-sorted as follows
B_y={B_y ^C,{tilde over (B)}_y ^LA} (3.98)
where B_y ^Cis the largest set of capacitive tree-branches, {tilde over (B)}_y ^LAis a temporary set of non-capacitive tree-branches.
Step 2: All link-branches in B_xcorresponding to trees in {tilde over (B)}_y ^LAare found. This set of link-branches is denoted as {tilde over (B)}_x ^LA. The combined set
{tilde over (B)}^LA={{tilde over (B)}_y ^LA,{tilde over (B)}_x ^LA} (3.99)
has no capacitive branches, moreover (3.99) can form a closed graph {tilde over (G)}=(Ñ,{tilde over (B)}^LA).
The remaining link-branches become
B _x ^LC =B _x −{tilde over (B)} _x ^LA (3.100)
Among the branches of B_x ^LC, there may be some which are inductive. These inductive link-branches can be identified and placed into set denoted B_x ^L ¹. After that, B_x ^LCcan be re-sorted as
B_x ^LC={B_x ^L ¹,B_x ^C} (3.101)
where B_x ^Cincludes the remaining non-inductive link-branches that may or may not be capacitive.

Step 3: The Minimum Spanning Tree Algorithm is Called

MinSTA({tilde over (G)},w _L)=
{tilde over (G)} _trees=(N,B _y ^LA) (3.102)
The corresponding set of link-branches is
B _x ^LA ={tilde over (B)} ^LA −B _y ^LA (3.103)
In the set B_x ^LA, all inductive branches are placed in the set of inductive link branches denoted as B_x ^L ². The remaining non-inductive links are
B _x ^A =B _x ^LA −B _x ^L ² (3.104)
Thereafter, the set B_x ^LAcan be reordered as
B_x ^LA={B_x ^A,B_x ^L ²} (3.105)
Step 4: B_x ^L ²is removed from the branch set {tilde over (B)}^LA. The remaining closed graph becomes Ĝ=({circumflex over (N)},{circumflex over (B)}^A). In order to accomplish this, it is necessary to identify all tree-branches in B_y ^LAthat are linked by B_x ^A. This branch set is denoted as B_y ^A. Thereafter, B_y ^LAcan be partitioned as
B_y ^LA{B_y ^L,B_y ^A} (3.106)
where B_y ^Lis found as usual
B _y ^L =B _y ^LA −B _y ^A (3.107)
It is important to note that none of the branches in {B_y ^A, B_x ^A} are inductive or capacitive.
Finally, the global branch set B can be re-ordered in terms of the smaller sets. Let the order be given as follows
B={B_y ^C,B_y ^L,B_y ^A,B_x ^A,B_x ^L ²,B_x ^L ¹,B_x ^C} (3.108)
The elementary networks are formed based on the branch sets as follows
{B_y ^C,B_x ^C}∈N_C (3.109)
{B_y ^L,B_x ^L ²,B_x ^L ¹}∈N_L (3.110)
{B_y ^A,B_x ^A}∈N_A (3.111)
Based on the branch order established in (3.108), the TCF of the node incidence matrix has the following structure
$\begin{matrix} A_{TCF} = [\begin{matrix} I_{τ \times τ}^{C} & 0 & 0 & {\hat{D}}_{τ \times t}^{CA} & {\hat{D}}_{τ \times h}^{L_{2}} & {\hat{A}}_{τ \times k}^{L_{1}} & {\hat{A}}_{τ \times z}^{C} \\ 0 & I_{η \times η}^{L} & 0 & 0 & {\hat{A}}_{η \times h}^{L} & 0 & 0 \\ 0 & 0 & I_{ζ \times ζ}^{A} & {\hat{A}}_{ζ \times t}^{A} & {\hat{C}}_{ζ \times h}^{LA} & 0 & 0 \end{matrix}] & (3.112) \end{matrix}$
where blocks have dimensions corresponding to the respective branch sets in (3.108). In particular

- τ=size (B_y ^C), as determined in (3.98)=
- η=size (B_y ^L), as determined in (3.107)
- ζ=size (B_y ^A), as determined in (3.106)=
- t=size (B_x ^A), as determined in (3.104)
- h=size(B_x ^L ²), as determined in (3.105)
- k=size (B_x ^L ¹), as determined in (3.101)
- z=size(B_x ^L), as determined in (3.101)

The previous procedures are the techniques for sequential identification of elementary networks. Depending on the relative order in which the networks are identified from the global graph, the sequence of steps in the procedure may differ. Similar procedures may be constructed in which the order of elementary network identification is different from the two cases considered above. However, as it is expected, the results of such procedures are topologically equivalent.

General Interconnected Networks

In the previous chapter, electrical circuits constructed of branches with components selected from one of the parameter sets (2.10)-(2.13) were considered. The topological conditions under which a circuit of a particular kind can be modeled as a single network of an appropriate type were also presented. Such differentiation was motivated on a basis of the type of equation and state variables (if any) required to describe the respective network. Furthermore, it was shown that as more freedom is allowed in terms of the circuit components and their topology, it is no longer possible to describe the entire circuit as a single type of network. Instead, for the circuits composed of branch models FIG. 4A, it was necessary to partition the global network N as
N=N_L∪N_A (3.113)
In the case of circuits constructed from branch models FIG. 1.2( b), the necessary partitioning was of the form
N=N_C∪N_A (3.114)
Finally, when both types of branch models FIGS. 4( a) and (b) were considered at the same time, it was shown that the circuit can be consistently described by viewing the corresponding global network as being partitioned as
N=(N _L ∪N _LC ∪N _C)∪N _A (3.115)
Even though (3.115) comes directly from the TCF (3.48), it may also be written in simplified form as
N=N_L∪N_C∪N_A (3.116)
Several observations can also be made about the actual coupling between DEs (3.72)-(3.73). The DE (3.72) has currents as state variables, and voltages as its driving input. The states in DE (3.73) are voltages, and it is the currents that represent the external driving force. Therefore, these two state equations are said to be input-output compatible. Extending this framework, models N_Land N_Ccan be viewed as interconnected dynamical systems each of which is represented by its own state equation that are coupled through their respective inputs and outputs and, for certain topologies, through an additional interconnection term due to the AE corresponding to N_A.

Generic Network Identification

Sometimes it may also be useful to partition the global system into several networks based not only on branch contents but also on some other additional constraints or criteria that may come from the description of the actual system and the topological layout of its sections. The property that must be maintained in any network partitioning is the ability to automatically produce KCL and KVL matrices that can be used to describe the corresponding networks and the coupling among them.
In constructing the KCL and KVL matrices, an algorithm for assembling the TCF is needed. Previously, in order to assemble the TCF, the multi-graph G_g=(N, B) was searched for a forest of spanning trees G_trees ^g=(N,B_y) with some particular topological properties reflected in an appropriate weight function. This very technique may be extended to accommodate other constraints in addition to those which come from branch components. For instance, it may be possible to have several, possibly pre-specified or guessed, networks. Then, a generic weight function associated with the branch set B can be defined as follows
$\begin{matrix} \begin{matrix} for \begin{matrix} \begin{matrix} j = 1 \dots q \\ ζ = L, C, A \end{matrix} \\ ξ = 1 \dots size (P_{ζ}) \end{matrix} & ω_{ζ}^{ξ} (b_{j}) = {\begin{matrix} 0, & if P_{ζ} {ξ} (j) = 0 \\ 1, & if P_{ζ} {ξ} (j) \neq 0 \\ \pm \infty, & if b_{j} \notin N_{ζ}^{k} \end{matrix} \end{matrix} & (3.119) \end{matrix}$
where P_ζ{ξ} is to be understood as ξ-th member of the corresponding parameter set P_ζ as defined in (2.10)-(2.13). If the member of P_ζ is a matrix, then by convention P_ζ{ξ} (j) is the j-th diagonal entry of the corresponding matrix. Similarly, when P_ζ{ξ} is a vector, then the notation P_{ζ{ξ} (j) is the j-th element. The N} _ζ ^kis the k-th network to be identified. The ± sign in front of the infinity symbol in (3.119) is taken such that the branch b_j∉N_ζ ^kis not likely to be included in G_trees ^g=(N,B_y) whether the search is performed for the MinST or MaxST.
Furthermore, it is desirable to assemble separate systems of DAEs along with the corresponding coupling terms for some, possibly pre-specified sections of a larger network. Ideally, the network should be partitioned in such a way that the number of variables coupling the corresponding systems of DAEs is minimized. If such a goal is feasible and the number of state variables in each of the DEs is significantly larger than the number of coupling variables, the corresponding networks may be viewed as being weakly coupled. In terms of the topology of weakly coupled networks, it is reasonable to expect that the number of common or connecting branches is small.
Thus, a generic step of partitioning the network N into two smaller ones may be considered as
N=N₁∪N₂ (3.120)
In doing so, the columns in the node incidence matrix would also be reorganized such that the corresponding TCF would have a particular block structure needed to assemble the KCL and KVL matrices on a network-by-network basis. In this regard, two types of structures of the right hand side of TCF have been heretofore encountered: lower-block triangular as in (3.8); and upper-block triangular as in (3.29). Even though these two TCFs were computed based on different topological objectives, they are in fact equivalent. Therefore, without loss of generality, for the network partitioning (3.120), the lower-block triangular TCF may be considered. Based of such structure for the TCF, the KCL and KVL matrices would have the following form
$\begin{matrix} A_{a} = A_{TCF} = [\begin{matrix} I_{1} & 0 & 0 & A_{1} \\ 0 & I_{2} & A_{2} & A_{21} \end{matrix}] & (3.121) \\ B_{b} = [\begin{matrix} 0 & - {(A_{2})}^{T} & I_{2} & 0 \\ - {(A_{1})}^{T} & - {(A_{21})}^{T} & 0 & I_{1} \end{matrix}] & (3.122) \end{matrix}$
If the partitioning (3.120) is final, then (3.121)-(3.122) would already require a particular choice of state variables for each network. That is, if a single set of state variables is to be used to completely describe a network, it is necessary to have a KCL or a KVL (whichever is appropriate for the selected states) to include only the branches corresponding to this network. With respect to (3.120), it appears from (3.121) that it is possible to write KCL including only the branches of N_J. Therefore, one can select the set of branches corresponding to the block Al in (3.121) to represent states such as independent currents. Such a choice of state variables would require N, to be an inductive network. Relative to network N₂, it is possible to write an independent KVL based on (3.122). This suggests a set of independent voltages to be the states and the network to be capacitive. If the TCF with its right-hand side being block diagonal (instead of block triangular) then the networks corresponding to these blocks would be topologically disconnected. In fact, for the multi-graph G^g=(N,B), the TCF would have as many diagonal blocks as there are graphs in G^g.
On the other hand, if partitioning (3.120) is not final, it is possible to recursively apply similar sub-network identification procedures to either or both N₁and N₂each time updating the TCF. Ideally, this procedure could be continued until a final collection of elementary networks is obtained. The result of this sequential network partitioning may be symbolically expressed as
N=N N . . . =N . . . u (3.123)
The corresponding TCF can be expressed in the general lower block-triangular
$\begin{matrix} A_{TCF} = [\begin{matrix} I_{1} & ⋮ & A_{1} \\ I_{2} & ⋮ & A_{2} & A_{21} \\ ⋱ & ⋮ & ⋰ & ⋮ & ⋮ \\ I_{k} & ⋮ & A_{k} & \dots & A_{k 2} & A_{k 1} \\ ⋱ & ⋮ & ⋰ & ⋮ & ⋮ & ⋮ \\ I_{n} & ⋮ & A_{n} & \dots & A_{nk} & \dots & A_{n 2} & A_{n 1} \end{matrix}] & (3.124) \end{matrix}$
Based on the TCF (3.124), and following the usual procedure, it is possible to assemble the KCL and KVL matrices for the global network as well as for each of the networks present in (3.123). Thus, based on (3.124), the KCL matrix for the k-th network (that is, the self-KCL matrix) is defined as
a [k Ak] (3.125)
Then, the corresponding KCL coupling matrix relating the k-th and i-th networks can be written from the k-th row of the TCF as
A_a ^ki=[0A_ki] (3.126)
Using definitions (3.125)-(3.126), the KCL equations for all networks in (3.123) can be expressed as
$\begin{matrix} A_{a}^{k} i_{br}^{k} + \sum_{i = 1}^{k - 1} A_{a}^{ki} i_{br}^{i} = 0, for k = 1, 2, \dots, n & (3.127) \end{matrix}$
where for k=1 the summation is defined to be zero. Furthermore, based on the structure of (3.126), it can be seen that only currents of the link branches of the other networks are needed in (3.127). Thus, defining new matrices as
D_a ^ki=A_ki (3.128)
the KCL equations (3.127) can be rewritten as
$\begin{matrix} A_{a}^{k} i_{br}^{k} + \sum_{i = 1}^{k - 1} D_{a}^{ki} i_{x}^{i} = 0, for k = 1, 2, \dots, n & (3.129) \end{matrix}$
The KVL matrices are defined in a similar way. In particular, the KVL self-matrix for the k-th network is determined from (3.124) to be
B _b ^k=[−(A _k)^T I _k] (3.130)
The KVL matrix coupling the k-th and i-th networks is assembled as
B _b ^ki=[−(A _ik)^T0 ] (3.131)
Thereafter, the KVL equations relating all networks in (3.123) can be written as
$\begin{matrix} B_{b}^{k} v_{br}^{k} + \sum_{i = k + 1}^{n} B_{b}^{ki} v_{br}^{i} = 0, for k = 1, 2, \dots, n & (3.132) \end{matrix}$
with the convention that for k=n, the summation is defined to be zero. Noting that only voltages corresponding to the tree branches of the other networks are needed in (3.132), the following definition can be made
C _b ^ki=−(A _ik)^T (3.133)
Then (3.132) can be rewritten with only the tree voltages under summation as
$\begin{matrix} B_{b}^{k} v_{br}^{k} + \sum_{i = k + 1}^{n} C_{b}^{ki} v_{y}^{i} = 0, for k = 1, 2, \dots, n & (3.134) \end{matrix}$
It can also be noted that (3.128) and (3.133) are related very much like the similar matrices in (3.71).
Furthermore, based on the KCL equations (3.129) and KVL equations (3.134) it is possible to express the state equations for general interconnected networks (3.123). In order to do this it is necessary to assume that there is no mutual coupling among branches belonging to different networks. As usual, the derivation goes through the dimensionality reduction procedure applied to each network in (3.123) that has state variables. Since there are two types of such networks, namely inductive and capacitive, there are only two types of state equations involved. In particular, if for some k, the network type index is ζ=L, meaning that network N_ζ ^k, is inductive, then the corresponding state equation is
$\begin{matrix} L_{x}^{k} \frac{\partial i_{x}^{k}}{\partial t} = - (R_{x}^{k} + \frac{\partial L_{x}^{k}}{\partial t}) i_{x}^{k} - B_{b}^{k} e_{br}^{k} - \sum_{i = k + 1}^{n} C_{b}^{ki} v_{y}^{i} & (3.135) \end{matrix}$
where all the reduced matrices are defined as before in terms of B_b ^k. Of course, in order to solve all DEs for multiple inductive networks simultaneously and without additional constraints, all state equations (3.135) must be input-output compatible. That is, for all i in (3.135) for which the network type index is ζ=L, it is necessary to have B_b ^ki=0, implying that the external forcing voltage may come only from other non-inductive networks. If it is not possible to satisfy this condition for two or more inductive networks in (3.123), then those networks should be combined into one. Similarly if, for some k, ζ=C, then network N_ζ ^kis capacitive, and the corresponding state equation is
$\begin{matrix} C_{y}^{k} \frac{\partial v_{y}^{k}}{\partial t} = - (G_{y}^{k} + \frac{\partial C_{y}^{k}}{\partial t}) v_{y}^{k} + A_{a}^{k} j_{br}^{k} - \sum_{i = 1}^{k - 1} D_{a}^{ki} i_{x}^{i} & (3.136) \end{matrix}$
Here, all reduced matrices are defined in terms of A_a ^k. Also, in order to be able to simultaneously model multiple capacitive networks without additional constraints, all state equations (3.136) must be input-output compatible. That is, all external forcing currents should come only from other non-capacitive networks, which implies that for all i in (3.135) for which ζ=C, it is necessary to have A_a ^ki=0. Whenever it is not possible to satisfy this condition for some capacitive networks in (3.123), then those networks should be re-combined. The relation of multiple elementary networks is symbolically depicted in FIG. 12. There, the external voltage and current sources are denoted as e_extand j_ext. Also, none of the inductive networks are connected with each other directly. Similarly, the capacitive networks are not coupled among each other through their branch variables. However, all networks are interconnected in such a way that each network with state variables has its proper inputs and outputs. On the other hand, algebraic networks may accommodate any inputs and outputs that are required by respective inductive and capacitive networks. Such an arrangement of networks allows all DAEs to be solved simultaneously.
The TCF (3.124) has a very general structure applicable to interconnected networks. The KCLs (3.127) and KVLs (3.132) also reflect possible coupling among all networks in (3.123) through their branch currents and voltages. If many of the networks in (3.123) are mutually de-coupled or weakly coupled, it can be expected that many of the block matrices with double subscripts in (3.124) are zero. For instance, if it is possible to make the right side of the TCF (3.124) block-bi-diagonal, the corresponding networks in (3.123) would be sequentially connected. Another useful way of formulating the relations among the networks is to have one (or maybe several) specific networks that represent all interconnections. In terms of the right side of the TCF, an attempt would be made to form a block diagonal structure as far down as possible by appropriately selecting the networks.

Selective Network Partitioning

When modeling a switched electrical network N=(G, P, S), it may be of practical interest to partition N into more than just three elementary networks N_L, N_C, and N_A. In general, the network partitioning may be of a more advanced form such as (3.123) where the branches may be assigned to networks based on some additional constraints. For instance, it is reasonable to establish additional constraints such that some of the networks of (3.123) remain unchanged throughout the entire simulation study. The same goal can also be pursued on the basis of elementary networks.
For convenience, it is assumed that each elementary network is partitioned into two smaller ones as
N_ζ=N_ζ ¹∪N_ζ ² (3.143)
where the running subscript is given as ζ=L, C, A. The network identification in (3.143) could be performed with the objectives given in the following subsections.
Identification of Networks with Variable Parameters
Some of the branch parameters of N_ζ may depend on time as well as on applied currents and voltages. In this case, the equations corresponding to N_ζ become nonlinear with time-varying coefficients. Therefore, it is desirable to partition the network as in (3.143) such that one of the networks, say N_ζ ¹, takes all branches with nonlinear and/or time varying parameters. It may be necessary to include other branches in order to make N_ζ ¹a proper network. After performing such network partitioning, it may become possible to assemble two systems of equations for N_ζ ¹and N_ζ ², respectively, such that each of them posses smaller dimensions. Then, the equations corresponding to N_ζ ²can be assembled once per topology, and the nonlinear equations with time-varying parameters corresponding to N_ζ ¹would be easier to deal with due to their reduction in size. In order to perform partitioning of N_ζ into N_ζ ¹and N_ζ ²such that their interconnection can be conveniently handled using dimensionality reduction procedures, it is necessary only to assign appropriate weights to the network branches as in (3.119).

Identification of Switched Networks

For large switched electrical networks, the automated procedure of equation generation becomes quite expensive. As a result of this computational expense, the models of systems with fast switching, such as models that include PWM, hysteresis current inverters, converters, etc., tend to require significantly longer CPU time to run the simulation. However, for a given network N=(G, P, S), it is likely that not all of its structural sections participate in the switching transients throughout the history of topological changes. One might, therefore, attempt to identify system sections with constant topology and group them into separate networks. This idea can also be applied to each elementary network of N. Therefore, with respect to (3.143), an elementary network N_ζ can be partitioned such that either N_ζ ¹or N_ζ ²posses a system of equations that need not be reassembled for each new topology. Thereafter, an attempt is made to exclude all such networks from the equation assembling procedures that are performed at each switching instance. Separating the global network N into its switching and non-switching parts would not only reduce the total amount of computations required per change in topology, but also provide a means for the local averaging of state equations for the switched subnetwork as is known in the art.

Network Partitioning for Stability Estimation

The structure of the DAEs produced by ASMG may also be utilized for the non-impedance-based stability analysis of energy conversion systems represented by their equivalent circuits. That is, instead of relying on the linearization of the state equations and using Nyquist-type criteria, it is possible to generate the DAEs in a form suitable for the Lyapunov analysis. For instance, if needed, a change of variables could be used to rewrite the state equations in the following autonomous form.
$\begin{matrix} \frac{\partial x_{L}}{\partial t} = f_{L} (x_{L}) + g_{L} (x_{L}, x_{C}) & (3.144) \end{matrix}$
where g_L(x_L, x_C) and g_C(x_L, x_C) are the interconnection terms that also include the algebraic network. The stability of (3.144)-(3.145) can be analyzed as similar to a problem of perturbed motion and studied using Lyapunov functions. However, for this technique, it is necessary to establish Lyapunov functions for each system in (3.144)-(3.145) ignoring the interconnections to begin with.
Alternatively, since it may be feasible to assemble a separate state equation for the part of the network N which is time invariant and linear, it may be advantageous to formulate the problem of stability analysis into a robust control framework. For instance, if it is possible to collect all nonlinear elements into a separate algebraic network with no current and/or voltage sources, then the terms that are due to N_Acan be viewed as a perturbation to the rest of the system.
The underlying power of the robust control framework is that it not only enables analysis of the stability of nonlinear systems in certain regions of the state space (attractors, invariant subspaces, etc.) but it also makes possible the design of controllers that ensure some desirable properties. In particular, considering quadratic Lyapunov functions, many analysis and design control problems associated with the model of N can be formulated as linear matrix inequalities (LMIs) that can be solved numerically. In this regard, it can be very advantageous to automatically assemble state and output equations for the global network N directly in a form permitting a linear differential inclusion (LDI) of a particular type.

Realization of the Global Network

In this section, a global network N=(G, P) with augmented topology is considered further. As shown in the previous chapter, a network of this kind, in general, possesses an interconnection of the form N=N_L∪N_C∪N_A. Such a network representation allows the corresponding systems of DAEs (3.72)-(3.75) to be assembled in a consistent manner for each elementary network with respective interconnections. For compactness, (3.72)-(3.75) can be written as
$\begin{matrix} \frac{\partial x_{LC}}{\partial t} = f (x_{LC}, y_{i υ}^{A}, u_{ej}, t) & (4.1) \\ g (x_{LC}, y_{i υ}^{A}, u_{ej}, t) = 0 & (4.2) \end{matrix}$
where the global state vector x_LC=[x_L, x_C] consists of minimal state vectors corresponding to the inductive and capacitive networks, y_iv ^A=[i_br ^A,v_br ^A] is a vector of branch currents and voltages corresponding to the algebraic network, and u_ej=[e_br,j_br] is an input vector that is composed of input voltage and current sources.
System (4.1) is a first-order system of ODEs, whereas (4.2) is a system of AEs. The AE (4.2) has the form (3.23), (3.25), (3.39), (3.41) and (3.74)-(3.75). If the elements of N_Aare some known, possibly nonlinear functions of applied currents and voltages, then (4.2) may not have a unique solution. In general, (4.2) may be solved using iterative methods for solving systems of nonlinear equations.
Global Network with Linear Parameters
Here, it is assumed that the network parameters are possibly time-varying, but do not depend on branch currents or voltages. Then, assuming that (3.75) is well-posed, the linear AE (4.2) always has a unique solution. Moreover, it is possible to solve (4.2) explicitly for the y_iv ^A. A After doing so, (4.1) and (4.2) can be merged into one system of first order ODEs. Thereafter, a system of DAEs describing the global network N can be mapped into a state space realization in its standard form
$\begin{matrix} \frac{\partial x_{LC}}{\partial t} = A (t) x_{LC} + B (t) u_{ej} & (4.3) \\ y_{i υ} = C (t) x_{LC} + D (t) u_{ej} & (4.4) \end{matrix}$
where the output vector is defined as y_iv=[i_br,v_br].
Before the state equations (3.72)-(3.73) can actually be written in the form of (4.3), some work with AEs (3.74)-(3.75) is needed. For simplicity of derivation, it is assumed that the state variables for the N_Land N_Care chosen to be independent currents and voltages as defined in (3.63) and (3.64), respectively. Then, defining matrix MA as
$\begin{matrix} M_{A} = [\begin{matrix} {\hat{A}}_{a}^{A} \\ {\hat{B}}_{b}^{A} {\hat{R}}_{br}^{A} \end{matrix}] & (4.5) \end{matrix}$
which is assumed to have full rank, (3.75) can be rewritten as
$\begin{matrix} \begin{matrix} {\hat{i}}_{br}^{A} = M_{A}^{- 1} [\begin{matrix} - C^{LA} \\ 0 \end{matrix}] i_{x} + M_{A}^{- 1} [\begin{matrix} 0 \\ - D^{CA} \end{matrix}] v_{y} + \\ M_{A}^{- 1} [\begin{matrix} 0 \\ - {\hat{B}}_{b}^{A} \end{matrix}] {\hat{e}}_{br}^{A} + M_{A}^{- 1} [\begin{matrix} 0 \\ - {\hat{B}}_{b}^{A} {\hat{R}}_{br}^{A} \end{matrix}] {\hat{j}}_{br}^{A} \\ = D_{I}^{A} i_{x} + D_{V}^{A} v_{y} + D_{e}^{A} {\hat{e}}_{br}^{A} + D_{j}^{A} {\hat{j}}_{br}^{A} \end{matrix} & (4.6) \end{matrix}$
The expression for the branch voltages (3.74) for N_Acan be rewritten as
$\begin{matrix} \begin{matrix} {\hat{v}}_{br}^{A} = {\hat{R}}_{br}^{A} M_{A}^{- 1} [\begin{matrix} - C^{LA} \\ 0 \end{matrix}] i_{x} + {\hat{R}}_{br}^{A} M_{A}^{- 1} [\begin{matrix} 0 \\ - D^{CA} \end{matrix}] v_{y} + \\ (I + {\hat{R}}_{br}^{A} M_{A}^{- 1} [\begin{matrix} 0 \\ - {\hat{B}}_{b}^{A} \end{matrix}]) {\hat{e}}_{br}^{A} + {\hat{R}}_{br}^{A} (I + M_{A}^{- 1} [\begin{matrix} 0 \\ - {\hat{B}}_{b}^{A} {\hat{R}}_{br}^{A} \end{matrix}]) {\hat{j}}_{br}^{A} \\ = C_{I}^{A} i_{x} + C_{V}^{A} v_{y} + C_{e}^{A} {\hat{e}}_{br}^{A} + C_{J}^{A} {\hat{j}}_{br}^{A} \end{matrix} & (4.7) \end{matrix}$
Substituting (4.6) into (3.72), and (4.7) into (3.73), and collecting terms, the state equations for N_Land N_Cbecome
$\begin{matrix} \begin{matrix} \frac{\partial i_{x}}{\partial t} = - L_{x}^{- 1} (R_{x} + \frac{\partial L_{x}}{\partial t} + C_{b}^{LA} C_{I}^{A}) i_{x} - L_{x}^{- 1} (C^{LC} + C_{b}^{LA} C_{V}^{A}) v_{y} - \\ L_{x}^{- 1} (B_{b}^{L} + C_{b}^{LA} C_{e}^{A}) e_{br}^{LA} - L_{x}^{- 1} C_{b}^{LA} C_{J}^{A} j_{br}^{A} \\ = A_{L}^{LA} i_{x} + A_{L}^{CA} v_{y} + B_{L}^{LA} e_{br}^{LA} + B_{L}^{A} j_{br}^{A} \end{matrix} & (4.8) \\ \begin{matrix} \frac{\partial v_{y}}{\partial t} = - C_{y}^{- 1} (G_{y} + \frac{\partial C_{y}}{\partial t} - D_{a}^{CA} D_{υ}^{A}) v_{y} - \\ C_{y}^{- 1} (D^{LC} + D_{a}^{CA} D_{I}^{A}) i_{x} + C_{y}^{- 1} (A_{a}^{C} + D_{a}^{CA} D_{J}^{A}) \\ j_{br}^{CA} - C_{y}^{- 1} D_{a}^{CA} D_{e}^{A} e_{br}^{A} \\ = A_{C}^{CA} v_{y} + A_{C}^{LA} i_{x} + B_{C}^{CA} j_{br}^{CA} + B_{C}^{A} e_{br}^{A} \end{matrix} & (4.9) \end{matrix}$
Finally (4.8) and (4.9) can be assembled to form a state equation in its standard form
$\begin{matrix} \frac{\partial}{\partial t} [\begin{matrix} i_{x} \\ v_{y} \end{matrix}] = [\begin{matrix} A_{L}^{LA} & A_{L}^{CA} \\ A_{C}^{LA} & A_{C}^{CA} \end{matrix}] [\begin{matrix} i_{x} \\ v_{y} \end{matrix}] + [\begin{matrix} B_{L}^{LA} & B_{L}^{A} \\ B_{C}^{A} & B_{C}^{CA} \end{matrix}] [\begin{matrix} e_{br} \\ j_{br} \end{matrix}] & (4.10) \end{matrix}$
The stage is set for assembling (4.4). Expressions (4.6) and (4.7) already provide currents and voltages for the branches of N_A. Thus, it is necessary to express currents and voltages for the branches of the two remaining networks, N_Land N_C. Starting with N_La vector of corresponding branch currents is determined from independent currents as
i _br ^L=(B _b ^L)^T i _x (4.11)
A vector of branch voltages for the N_Lcan be written as
$\begin{matrix} \begin{matrix} v_{br}^{L} = R_{br} i_{br}^{L} + L_{br} \frac{\partial i_{br}^{L}}{\partial t} + \frac{\partial L_{br}}{\partial t} i_{br}^{L} + e_{br}^{L} \\ = (R_{br} + \frac{\partial L_{br}}{\partial t}) i_{br}^{L} + L_{br} \frac{\partial i_{br}^{L}}{\partial t} + e_{br}^{L} \end{matrix} & (4.12) \end{matrix}$
Substituting (4.11) into (4.12), an expression for branch voltages in terms of independent currents can be obtain as
$\begin{matrix} v_{br}^{L} = (R_{br} + \frac{\partial L_{br}}{\partial t}) {(B_{b}^{L})}^{T} i_{x} + {L_{br} (B_{b}^{L})}^{T} \frac{\partial i_{x}}{\partial t} + e_{br}^{L} & (4.13) \end{matrix}$
After substituting state equation (4.8) into (4.13) and collecting terms, the final expression for branch voltages becomes
$\begin{matrix} v_{br}^{L} = ((R_{br} + \frac{\partial L_{br}}{\partial t}) {(B_{b}^{L})}^{T} - {L_{br} (B_{b}^{L})}^{T} L_{x}^{- 1} (R_{x} + \frac{\partial L_{x}}{\partial t} + C_{b}^{LA} C_{I}^{A})) i_{x} - ({L_{br} (B_{b}^{L})}^{T} L_{x}^{- 1} (C^{LC} + C_{b}^{LA} C_{V}^{A})) v_{y} - ({L_{br} (B_{b}^{L})}^{T} L_{x}^{- 1} (B_{b}^{L} + C_{b}^{LA} C_{V}^{A}) - I^{L}) e_{br}^{LA} - ({L_{br} (B_{b}^{L})}^{T} L_{x}^{- 1} C_{b}^{LA} C_{j}^{A}) j_{br}^{A} = C_{L}^{LA} i_{x} + C_{L}^{CA} v_{y} + D_{L}^{LA} e_{br} + D_{L}^{A} j_{br} & (4.14) \end{matrix}$
where I^Lis an identity-like matrix with ones in diagonal entries corresponding to branches in N_Land zeros elsewhere.
The equations for a, capacitive network are assembled in a similar way. That is, recalling (3.64), a vector of corresponding branch voltages is determined as
v _br ^C(A _a ^C)^T v _y (4.15)
A vector of branch currents for the N_Ccan be found as
$\begin{matrix} i_{br}^{C} = (G_{br} + \frac{\partial C_{br}}{\partial t}) v_{br}^{C} + C_{br} \frac{\partial v_{br}^{C}}{\partial t} - j_{br}^{C} & (4.16) \end{matrix}$
Combining (4.15) and (4.16) yields
$\begin{matrix} i_{br}^{C} = (G_{br} + \frac{\partial C_{br}}{\partial t}) {(A_{a}^{C})}^{T} v_{y} + {C_{br} (A_{a}^{C})}^{T} \frac{\partial v_{y}}{\partial t} - j_{br}^{C} & (4.17) \end{matrix}$
Similar to (4.14), (4.9) can be substituted for the state derivative in (4.17). Collecting corresponding terms, the following result is obtained
$\begin{matrix} i_{br}^{C} = ((G_{br} + \frac{\partial C_{br}}{\partial t}) {(A_{a}^{C})}^{T} - {C_{br} (A_{a}^{C})}^{T} C_{y}^{- 1} (G_{y} + \frac{\partial C_{y}}{\partial t} - D_{a}^{CA} D_{V}^{A})) v_{y} + ({C_{br} (A_{a}^{C})}^{T} C_{y}^{- 1} (D^{LC} + D_{a}^{CA} D_{I}^{A})) i_{x} + ({C_{br} (A_{a}^{C})}^{T} C_{y}^{- 1} (A_{a}^{C} + D_{a}^{CA} D_{j}^{A}) - I^{C}) j_{br}^{CA} + ({C_{br} (A_{a}^{C})}^{T} C_{y}^{- 1} D_{a}^{CA} D_{e}^{A}) e_{br}^{A} = C_{C}^{CA} v_{y} + C_{C}^{LA} i_{x} + D_{C}^{CA} j_{br}^{CA} + D_{C}^{A} e_{br}^{A} & (4.18) \end{matrix}$
where I^Cis also an identity-like matrix with ones in diagonal entries corresponding to branches in N_Cand zeros elsewhere.
Assuming that vectors of currents and voltages for all branches of the global network N can be assembled (concatenated) from the corresponding currents and voltages of N_L, N_C, and N_Aas
i _br =i _br ^L +i _br ^C +i _br ^A (4.19)
v _br =v _br ^L +v _br ^C +v _br ^A (4.20)
the expressions (4.14) and (4.18) are put together to form (4.4) which becomes
$\begin{matrix} [\begin{matrix} i_{br} \\ v_{br} \end{matrix}] = [\begin{matrix} {(B_{b}^{L})}^{T} + C_{C}^{LA} + D_{I}^{A} & C_{C}^{CA} + D_{V}^{A} \\ C_{C}^{LA} + C_{I}^{A} & {(A_{a}^{C})}^{T} + C_{L}^{CA} + C_{V}^{A} \end{matrix}] [\begin{matrix} i_{x} \\ v_{y} \end{matrix}] + [\begin{matrix} D_{C}^{A} + D_{e}^{A} & D_{C}^{CA} + D_{J}^{A} \\ D_{L}^{LA} + C_{e}^{A} & D_{L}^{A} + D_{j}^{A} \end{matrix}] [\begin{matrix} e_{br} \\ j_{br} \end{matrix}] & (4.21) \end{matrix}$
Thus, (4.10) and (4.21) form an automatically assembled minimal realization for the global network N.
Global Network with Nonlinear Parameters
Heretofore, only networks with linear time-varying parameters were considered. Electrical networks with such parameters are linear in the sense that it is possible to apply superposition with respect to state variables. Such networks can be successfully modeled using topological algorithms and dimensionality reduction procedures and that, if needed, the resulting system of DAEs can be mapped into a standard form state space realization. In general, network parameters may not only depend on time but also be some functions of state variables such as inductor currents and fluxes, and capacitor voltages and charges. For instance, a branch resistance may depend on the current flowing through or the voltage applied to the branch. Also, if nonlinear magnetic properties are to be included in the model, the corresponding branch inductances become state dependent. This phenomena is caused by the saturation of the magnetic materials used in inductors. In modeling semiconductor power electronic devices, the effective capacitance of the junction may also depend on voltage or current depending upon the type of device considered. Therefore, it should be considered that
R _br =R _br(t,x _L) (4.22)
G _br =G _br(t,x _C) (4.23)
L _br =L _br(t,x _L) (4.24)
C _br =C _br(t,x _C) (4.25)
As before, a system of DAEs for each elementary network and their interconnection will be derived ignoring inactive branches for the sake of convenience. First, the system of AEs (3.75) for the network N_Ais rewritten as
$\begin{matrix} [\begin{matrix} {\hat{A}}_{a}^{A} \\ {\hat{B}}_{b}^{A} {\hat{R}}_{br}^{A} \end{matrix}] {\hat{i}}_{br}^{A} - [\begin{matrix} - C^{LA} & 0 \\ 0 & - D^{CA} \end{matrix}] [\begin{matrix} i_{x} \\ v_{y} \end{matrix}] - [\begin{matrix} 0 & 0 \\ - {\hat{B}}_{b}^{A} & - {\hat{B}}_{b}^{A} {\hat{R}}_{br}^{A} \end{matrix}] [\begin{matrix} {\hat{e}}_{br}^{A} \\ {\hat{j}}_{br}^{A} \end{matrix}] = 0 & (4.26) \end{matrix}$
which is a system of nonlinear equations that can be solved iteratively for î_br. After doing so, the vector of corresponding branch voltages, {circumflex over (v)}_bris computed using (3.74).
To assemble DAEs for the inductive network, the voltage equation for branches of this network may be expressed in its most general form
$\begin{matrix} v_{br}^{L} = R_{br}^{L} i_{br}^{L} + \frac{\partial λ_{br}}{\partial t} + e_{br}^{L} & (4.27) \end{matrix}$
Taking into account other networks KVL (3.70) can be written as
B _b ^L v _br ^L+Σ^CA=0 (4.28)
where the symbol Σ^CArepresents all terms that are due to interconnection of N_Lwith capacitive and algebraic networks.
Working with nonlinear inductances, (4.24) represents an inconvenience in that it becomes necessary to compute the total time derivative of the matrix L_xwhich, in turn, contains partial derivatives of L_xwith respect to the vector i_x. A partial derivative of a matrix with respect to vector is a “large” and cumbersome object. However, the nonlinear relationship between fluxes and currents must be incorporated into the system of DAEs for the network N_L. Therefore, instead of (4.24), the relationship between fluxes and currents can be modeled using a saturation function that would “penalize” the flux linkages at high currents. Using this technique of representing magnetic saturation, the branch fluxes are expressed as
λ_br =L _br i _br ^L−φ _br(i _br ^L) (4.29)
where φ_br, is a vector-valued saturation function. From this point, (4.29) is substituted into (4.27) and the dimensionality reduction procedure using KVL (4.28) is applied to the resulting voltage equation. In doing so, the vector of independent fluxes is found as follows
λ_x =L _x i _x −B _b ^Lφ_br[(B _b ^L)^T i _x ]=L _x i _x−φ_x(i _x) (4.30)
where, by analogy to all other reduced quantities, φ_x(i_x) is the reduced saturation function.
A direct result of the dimensionality reduction procedure is a state equation with independent fluxes as state variables. For consistency, this equation is written as
$\begin{matrix} \frac{\partial λ_{x}}{\partial t} = - R_{x} i_{x} - B_{b}^{L} e_{br}^{L} - Σ^{CA} & (4.31) \end{matrix}$
Equations (4.30)-(4.31) form a system of DAEs for the network N_L. Here, it is necessary to solve the nonlinear AE (4.30) for i_xat each integration step, or more precisely, at each call to the derivative function (4.31). In general, iterative methods would be used for solving (4.30) for the vector of currents. For small integration steps, Newton's method with an initial guess taken as i_xfrom the previous integration step could yield fast convergence. It may also be possible to re-parameterize the saturation function φ_x(i_x) into φ′_x(λ_x). If this is achieved, then (4.30) would become
L _x i _x=λ_x+φ′_x(λ_x) (4.32)
which is a linear equation with respect to i_x.
Instead of re-parameterization of the saturation function φ_x(i_x), i_xcan be chosen to be the state vector in order to avoid having to solve a system of nonlinear equations. Using the chain rule, the time derivative of fluxes in (4.30) can be expressed as
$\begin{matrix} \frac{\partial λ_{x}}{\partial t} = \frac{\partial λ_{x}}{\partial i_{x}} \frac{\partial i_{x}}{\partial t} + \frac{\partial λ_{x}}{\partial t} = (L_{x} - \frac{\partial ϕ_{x}}{\partial i_{x}}) \frac{\partial i_{x}}{\partial t} + \frac{\partial L_{x}}{\partial t} i_{x} & (4.33) \end{matrix}$
Combining (4.30) and (4.33), the state equation for the network N_Lwith currents as the state variables can be written in implicit form as
$\begin{matrix} (L_{x} - \frac{\partial ϕ_{x}}{\partial i_{x}}) \frac{\partial i_{x}}{\partial t} = - (R_{x} + \frac{\partial L_{x}}{\partial t}) i_{x} - B_{b}^{L} e_{br}^{L} - Σ^{CA} & (4.34) \end{matrix}$
The vector of branch currents is computed as usual using (4.11), and (4.28) can be used to compute branch voltages. Applying the chain rule to (4.30) and utilizing (4.11), the derivative of branch fluxes can be expressed as
$\begin{matrix} \frac{\partial λ_{br}}{\partial t} = (L_{br} - \frac{\partial ϕ_{br}}{\partial i_{br}}) {(B_{b}^{L})}^{T} (\frac{\partial i_{x}}{\partial t}) + \frac{\partial L_{br}}{\partial t} {(B_{b}^{L})}^{T} i_{x} & (4.35) \end{matrix}$
Combining (4.28) and (4.35), the output equation becomes
$\begin{matrix} v_{br}^{L} = (L_{br} - \frac{\partial ϕ_{br}}{\partial i_{br}}) {(B_{b}^{L})}^{T} (\frac{\partial i_{x}}{\partial t}) + (R_{br}^{L} + \frac{\partial L_{br}}{\partial t}) {(B_{b}^{L})}^{T} i_{x} + e_{br}^{L} & (4.36) \end{matrix}$
In (4.36), there is a term proportional to the derivative of the state vector. If (4.35) can be expressed in explicit form, it could be substituted into (4.36). However, since the derivative of the state vector is computed anyway, this term can be present in (4.36) without any computational disadvantage.
A similar system of DAEs can also be assembled for capacitive networks. The current equation for branches of N_Cis
$\begin{matrix} i_{br}^{C} = G_{br}^{C} v_{br}^{C} + \frac{\partial q_{br}}{\partial t} - j_{br}^{C} & (4.37) \end{matrix}$
and KCL (3.67) can be written here as
A _a ^C i _br ^C+Σ^LA=0 (4.38)
For the capacitive network, the nonlinear relation between capacitor charges and voltages can be represented in a form similar to magnetic saturation. That is
q _br =C _br v _br ^C−φ_br(v _br ^C) (4.39)
where φ_bris also a vector-valued function representing a part of capacitance that depends on the voltage. Using the dimensionality reduction procedure with (4.15) and (4.38), a vector of independent charges can be found as
q _y =C _y v _y −A _a ^Cφ_br((A _a ^C)^T v _y)=C _y v _y−φ_y(v _y) (4.40)
Using the chain rule and (4.15), the derivatives of (4.39)-(4.40) with respect to time are found as
$\begin{matrix} \frac{\partial q_{br}}{\partial t} = (C_{br} - \frac{\partial φ_{br}}{\partial v_{br}}) {(A_{a}^{C})}^{T} (\frac{\partial v_{y}}{\partial t}) + \frac{\partial C_{br}}{\partial t} {(A_{a}^{C})}^{T} v_{y} & (4.41) \\ \frac{\partial q_{y}}{\partial t} = (C_{y} - \frac{\partial φ_{y}}{\partial v_{y}}) \frac{\partial v_{y}}{\partial t} + \frac{\partial C_{y}}{\partial t} v_{y} & (4.42) \end{matrix}$
Finally, following the procedure, the state equation for N_Ccan also be obtained in implicit form as
$\begin{matrix} (C_{y} - \frac{\partial φ_{y}}{\partial v_{y}}) \frac{\partial v_{y}}{\partial t} = - (C_{y} + \frac{\partial C_{y}}{\partial t}) v_{y} + A_{a}^{C} j_{br}^{C} - Σ^{LA} & (4.43) \end{matrix}$
The vector of branch currents is computed based on (4.37) and (4.41) as
$\begin{matrix} i_{br}^{C} = (C_{br} - \frac{\partial φ_{br}}{\partial v_{br}}) {(A_{a}^{C})}^{T} (\frac{\partial v_{y}}{\partial t}) + (G_{br}^{C} + \frac{\partial C_{br}}{\partial t}) {(A_{a}^{C})}^{T} v_{y} - j_{br}^{C} & (4.44) \end{matrix}$
The equations developed thus far can be written more compactly. For the algebraic network, (4.26) and (3.74) have the following form
f _A(î _br ^A , i _x , v _y ,ê _br ^A ,ĵ _br ^A ,t)=0 (4.45)
{circumflex over (v)} _br ^A =q _A(î _br ^A ,i _x ,v _y ,ê _br ^A ,ĵ _br ^A ,t) (4.46)
The two state equations (4.34) and (4.43) are rewritten as
$\begin{matrix} M_{L} (t) \frac{\partial i_{x}}{\partial t} = f_{L} (i_{x}, v_{y}, v_{br}^{A}, e_{br}^{L}, t) & (4.47) \\ M_{C} (t) \frac{\partial v_{x}}{\partial t} = f_{C} (i_{x}, v_{y}, j_{br}^{A}, j_{br}^{C}, t) & (4.48) \end{matrix}$
where M_L(t) and M_C(t) are so-called mass matrices that can dependent on time and state.
Combined, (4.11) and (4.38) are written compactly as
$\begin{matrix} y_{br}^{L} = g_{L} (\frac{\partial i_{x}}{\partial t}, i_{x}, e_{br}^{L}, t) & (4.49) \end{matrix}$
Similarly, (4.44) together with (4.15) are expressed as
$\begin{matrix} y_{br}^{C} = g_{C} (\frac{\partial v_{y}}{\partial t}, v_{y}, j_{br}^{C}, t) & (4.50) \end{matrix}$
where y_br ^L=[i_br ^L,v_br] and y_br ^C=[i_br ^C,v_br ^C] are the output vectors.
Equations (4.45)-(4.50) form a system of DAEs which describe the global network N. Here, (4.45) is a system of nonlinear equations that could be solved numerically. Then, (4.46) and (4.49)-(4.50) are explicit systems of AEs that are not expensive to evaluate. The implicit systems of DEs (4.47)-(4.48) may be solved for the respective derivatives using some efficient techniques developed for systems of linear equations. In such arrangements, a linear solver would be called at each call to the corresponding derivative function. Alternatively, there are efficient numerical techniques that have been designed to accommodate time dependent mass matrices.

Switched Electrical Networks

A network of a general kind was previously defined as N=(G, P, s), where s is some given topology. If N is a switched network, the vector s is different for each new topology. In order for the constructed model to be meaningful, N should be a finite electrical network for each encountered topology. Throughout the simulation process, the network N may switch sequentially from an initial topology s₀to the final topology s, say as
s₀, s₁, . . . , s_i, . . . , s_j, . . . , s_m (5.1)
Based on this sequence, a so-called topology matrix can be formed as
S=[s₀, s₁, . . . , s_i, . . . , s_j, . . . , s_m] (5.2)
For each distinct topology of (5.1), there is a minimal realization for the inductive and capacitive networks present in N. In general, these minimal realizations may have a different number of states for each topology. That is, for a given topology s_iof the sequence (5.1), the state vectors of N_Land N_Cbelong to some appropriate vector spaces given topology s_iof the sequence (5.1), the state vectors of N_Land N_Cbelong to some appropriate vector spaces
x_L ⁱ∈
^α ⁱ (5.3)
x_C ⁱ∈
^β ⁱ (5.4)
where the dimensions of the minimal state vectors for the networks N_Land N_Care denoted with respect to (5.1) as
α₀, α₁, . . . , α_i, . . . , α_j, . . . , α_m (5.5)
β₀, β₁, . . . , β_i, . . . , β_j, . . . , β_m (5.6)
respectively.
It is also expected that in the topological sequence (5.1), there exists an s_iand S_jin which the number of state variables is maximum for the realizations of N_Land N_C, respectively. Of course, there may be several distinct topologies that require the same maximum number of states for their respective realizations. Assuming that it is always possible to obtain a non-minimal realization of a desired size by including some redundant states, it is possible to define the maximum numbers of state variables required to implement networks N_Land N_Cas
α=max{α₀, α₁, . . . , a_i, . . . , α_j, . . . , α_m} (5.7)
β=max{β₀, β₁, . . . , β_i, . . . , β_j, . . . , β_m} (5.8)
respectively.
Furthermore, since for each single topology s_ithe network N=(G, P, s_i) must be finite, the same is required for the switched network N=(G, P, S). Symbolically this requirement can be rewritten as
N=(G,P,S)∈N _q ^α+β (5.9)
where α and β are defined in (5.7) and (5.8), and α+β≦q. Sometimes α and β are referred to as the network complexity. Also, since S is a topology matrix, N=(G, P, S) is equivalent to a family of networks.
The sequence (5.1) may contain repeated topologies, and the actual number of distinct topologies may be small relative to the length m of the complete sequence of topologies throughout the entire simulation history. In this case, a reduced topology matrix that includes only the distinct topological vectors may be defined. Without loss of generality, it can be assumed that in the sequence s₀, s₁, . . . , s_i, . . . s_j, . . . s_r, . . . , s_mfor some r≦m, the first r+1 vectors are distinct. Thereafter, a reduced topology matrix is defined as
S_r=[s₀, s₁, . . . , s_r] (5.10)
With respect to the network N, it is possible to say that S_rspans all encountered allowable topologies.
Calculation of Initial Conditions after Switching Event
Here, the two types of switched networks, namely inductive and capacitive are considered. For this case, it is required that the currents through inductive branches and the voltages across capacitive branches are continuous as the network undergoes topological changes. Then, the continuity conditions for inductive and capacitive networks can be written as
i _i ^L =i _i+1 ^L =i _br ^L, and ∥i _br ^L∥_∞<∞ (5.11)
v _i ^C =v _i+1 ^C =v _br ^C, and ∥v _br ^C∥_∞<∞ (5.12)
In other words, the vectors (or more precisely trajectories) i_br ^Land v_br ^Cmust be bounded and continuous across topological boundaries. Recalling how i_br ^Land v_br ^Care related to the vectors of independent inductor currents and capacitor voltages, (5.11)-(5.12) can also be rewritten as
(B _i ^L)^T i _x ⁱ=(B _i+1 ^L)^T i _x ⁱ⁺¹ =i _br ^L (5.13)
(A _i ^C)^T v _y ⁱ=(A _i+1 ^C)^T v _y ⁱ⁺¹ =v _br ^C (5.14)
The KVL and KCL matrices corresponding to the second topology s_i+1can also be expressed in terms of the branch order in the respective TCF as
B _i+1 ^L=└−(Â _i+1 ^L)^T I _i+1 ^L0┘(T _L ⁱ⁺¹)^T (5.15)
A _i+1 ^L =I _i+1 ^C Â _i+1 ^L0┘(T _C ⁱ⁺¹)^T (5.16)
where T_L ⁱ⁺¹and T_C ⁱ⁺¹are appropriate permutation matrices. Based on the structure of (5.15)-(5.16), it is possible to define the following right pseudo inverses
(B _i+1 ^L)⁺ =T _L ⁱ⁺¹[0I _i+1 ^L0]^T=(B _i+1 ^base)^T (5.17)
(A _i+1 ^L)⁺ =T _C ⁱ⁺¹ [I _i+1 ^C00]^T=(A _i+1 ^base)^T (5.18)
It can be noted that B_i+1 ^baseand A_i+1 ^baseare full-rank matrices containing only columns that make the basis in either case. Therefore, (5.17)-(5.18) can be used to compute the corresponding initial conditions for the new topology. In particular, the initial values of the independent currents and voltages for the topology s_i+1, can be computed as
i _x ⁱ⁺¹ =B _i+1 ^base i _i ^L (5.19)
v _y ⁱ⁺¹ =A _i+1 ^base v _i ^C (5.20)
That is, based on branch currents and voltages just before commutation, it is possible to compute initial values for the state variables for the next topology such that the usual continuity conditions are satisfied.
In the present embodiment, the values in the system are constrained to represent finite currents and voltages. Thus, in resistive networks,
∥i _i ^A∥_∞<∞, and ∥i _i+1 ^A∥ _∞<∞ (5.24)
∥v _i ^A∥_∞<∞, and ∥v _i+1 ^A∥_∞<∞ (5.25)
In networks with state variables, this constraint may be expressed as
i _i ^L =i _i+1 ^L =i _br ^L, and ∥i _br ^L∥_∞<∞ (5.26)
v _i ^C =v _i+1 ^C =v _br ^C, and ∥v _br ^C∥_∞<∞ (5.27)

Non-Minimal Realization

Thus far, all state equations for the networks N_Land N_Cwere assembled algorithmically based on topological conditions using the minimal number of state variables required to represent the dynamics of each network. That is, for each individual topology, the system of DAEs constitutes a minimal state-space realization. Since the topology of N=(G, P, s) is changing, it is reasonable to expect that the dimensions of the minimal state vectors may also change. In general, the size of the minimal realization for each network N_Lor N_Cmay vary significantly throughout the sequence of topologies s₀, s₁, . . . , s_i, . . . , s_m. On the other hand, for analytical and/or numerical reasons, it may be convenient to have matrices and vectors that do not change their sizes dynamically from one topology to the next. For this reason, in a computer implementation of the ASMG, matrices A_a ^Cand B_b ^Lcan be defined so that the dimensions do not change throughout the entire simulation of N=(G, P, S). Whenever needed, the appropriate number of zero rows can be appended to the bottom in order to maintain the same matrix size. That is, A_a ^Cand B_b ^Lcan be modified as
$\begin{matrix} {\tilde{A}}_{a}^{C} = [\begin{matrix} A_{a}^{C} \\ 0 \end{matrix}] & (5.36) \\ {\tilde{B}}_{b}^{L} = [\begin{matrix} B_{b}^{L} \\ 0 \end{matrix}] & (5.37) \end{matrix}$
The total number of rows for A_a ^Cand for B_b ^Lis then determined by the maximum number of state variables (5.7)-(5.8) needed to represent networks N_Land N_Cfor any topology in the sequence s₀, s₁, . . . , s_i, . . . , s_m. Furthermore, Ã_a ^Cand for {tilde over (B)}_b ^Lusing in the dimensionality reduction procedures for N_Land N_C, respectively, the reduced matrices L_xand C_ywill be block-diagonal, with full-rank upper-left blocks and zeros elsewhere. Therefore, having such a structure, these matrices are also block-invertible. Moreover, the same block-diagonal structure is preserved under any non-singular coordinate transformation applied to the state variables in (2.106) or (2.111).
Besides having the same state-space dimensions for all encountered topologies, it is also possible to make some specific state variables redundant. For a network N_ζ with state variables and constant parameters, a minimal state-space realization can be assembled in its standard form
$\begin{matrix} \frac{\partial x_{ζ}}{\partial t} = A_{ζ} x_{ζ} + B_{ζ} u y_{iv} = C_{ζ} x_{ζ} + D_{ζ} u & (5.38) \end{matrix}$
where ζ=L, C, meaning that the network N_ζmay be inductive or capacitive.
A coordinate transformation matrix can be defined as
$\begin{matrix} K_{ζ} = [\begin{matrix} I \\ M \end{matrix}] & (5.39) \end{matrix}$
and the new vector of state variables as
{tilde over (x)}_ζ=K_ζx_ζ (5.40)
where M can be any real matrix of appropriate dimensions. Note that Mx_ζ is an additional vector of redundant states. That is, any element of this vector is nothing more than some linear combination of the state variables in x_ζ. Then, defining the right pseudo inverse of K_ζ to be of the form
K_ζ ^†=[I0] (5.41)
(5.38) can be rewritten in terms of the new states as
$\begin{matrix} \frac{\partial \tilde{x}}{\partial t} ζ = K_{ζ} A_{ζ} K_{ζ}^{†} {\tilde{x}}_{ζ} + K_{ζ} B_{ζ} u = {\tilde{A}}_{ζ} {\tilde{x}}_{ζ} + {\tilde{B}}_{ζ} u y_{iv} = C_{ζ} K_{ζ}^{†} {\tilde{x}}_{ζ} + D_{ζ} = {\tilde{C}}_{ζ} {\tilde{x}}_{ζ} + {\tilde{D}}_{ζ} u & (5.42) \end{matrix}$
This transformation of variables can be carried through it, instead of (5.38), the more general system of DAEs (4.10), (4.21), and (4.45)-(4.50) is considered.

Minimal Common Realization

Herein, a switched electrical network N=(G, P, S)∈N_q ^α+β with the sequence of topologies given in (5.1) is considered. It is assumed that N is switching between two distinct generic topologies, s_iand s_i+1. As usual, for any single topology s, the global network N=(G, P, s) possesses an interconnection of elementary networks of the form N=N_L∪N_C∪N_A. If inactive branches are included in their respective elementary networks, the commutation conditions (5.24)-(5.27) must also be satisfied with respect to some fixed branch order. In particular, for networks with state variables, the continuity of inductor currents and capacitor voltages across topological boundaries results in continuity of respective branch currents and voltages as expressed in (5.11)-(5.12).
In addition to the switching of the branch currents and voltages, it is also interesting to consider the commutation of state variables between topologies. For this purpose, a network with state variables N_ζ, where ζ=L, C, is considered. If N_ζis modeled using its respective minimal state equation, then it is reasonable to expect that the state vector x_ζwill also change its dimension from one topology to the next. Specifically, based on (5.3)-(5.4) and with respect to s_iand s_i+1, the state vector x_ζ can be related to the corresponding vector space as
x_ζ ⁱ∈
^γ ⁱ (5.43)
x^ζ ⁱ⁺¹∈
^γ ⁱ⁺¹ (5.44)
where the index γ_i=α_i, β_idenotes the size of the minimal state vector for inductive or capacitive network as given in (5.5)-(5.6) with respect to the sequence (5.1). Naturally, since the topologies s_iand s_i+1, are distinct, with regard to (5.43)-(5.44), it is reasonable to expect that
γ_i≠γ_i+1 (5.45)
and
x _ζ ⁱ ≠x _ζ ⁱ⁺¹ (5.46)
from which it follows that for the simulation process of switched network N_ζ, the state variables are discontinuous. Note that (5.46) is likely to be the case not only because of (5.45), but also due to the fact that the natural state variables (independent currents, fluxes, voltages, and charges) are selected algorithmically based only on the topological information for the active portion of the network. That is, the state vector x_ζ ⁱis assembled based only on N=(G, P, s_i). Similarly, the state variables in x_ζ ⁱ⁺¹are selected based on knowledge of N=(G, P, s_i+1). Therefore, even though the dimensions of x_ζ ⁱand x_ζ ⁱ⁺¹may happen to be the same, in general, the state variables will not be continuous across topological boundaries.
As shown in the preceding section, it is always possible to obtain a non-minimal realization for each topology of the sequence (5.1) such that instead of (5.43)-(5.44),
x_ζ ⁱ, x_ζ ⁱ⁺¹∈
(5.47)
where “γ” α, β denotes the maximum required size for the state vectors as defined in (5.7)-(5.8). Realizations of the maximum required size for which (5.47) holds can be obtained by appropriately applying the transformation of variables (5.39)-(5.41) for each new topology of the sequence (5.1). However, in general, condition (5.47) does not guarantee continuity of state variable across the topological boundaries, and {tilde over (x)}_ζ ⁱ∉C.
In addition to (5.47), it is desirable to replace (5.46) with an expression of the form
{tilde over (x)} _ζ ⁱ ={tilde over (x)} _ζ ⁱ⁺¹ (5.48)
If the transformation of state variables is performed in such a way that (5.39)-(5.41) holds for every switching event in the sequence (5.1), the resulting model would exhibit global state-space continuity, x_ζ∈C. The ASMG with the state scheduling algorithm developed in this thesis is capable of mapping each network incidence N=(G, P, s_i) into a system of DAEs. However, there are many numerical as well as analytical reasons that make the property x_ζ∈C very desirable.
Without condition (5.48), the simulation of the network N=(G, P, S) is essentially a concatenation of solutions of m+1 initial value problems (IVPs) corresponding to the time intervals between topologies (5.1). Therefore, the simulation would require re-computing the initial conditions and restarting the integration routine for each new topology. In the case of continuous state variables, the need for re-initializing the integration routine would disappear.
If conditions (5.47)-(5.48) are satisfied for all switching instances of (5.1), the corresponding systems of DAEs for different topologies also become compatible. That is, having state continuity x_ζ∈C, it is possible to apply averaging techniques known in the art to the state equations obtained for distinct topologies. Another important advantage of working with state-space realizations for which (5.47)-(5.48) are enforced, is that it becomes possible to apply Lyapunov-based stability analysis for polytopic systems. Therefore, the development of a deterministic state selection algorithm subject to the continuity constraints (5.47)-(5.48) should be considered.
The first step in deriving the state selection algorithm subject to {tilde over (x)}_ζ∈C is to define the change of variables (5.39)-(5.41) for two adjacent generic topologies s_iand s_i+1, and enforce (5.47)-(5.48). To achieve this, it is necessary to relate state vectors x_ζ ⁱand x_ζ ⁱ⁺¹through some matrices such that
x _ζ ⁱ=ⁱ T _ζ ⁱ⁺¹ x _ζ ⁱ⁺¹ (5.49)
x _ζ ⁱ⁺¹=ⁱ⁺¹ T _ζ ⁱ x _ζ ⁱ (5.50)
where ⁱT_ζ ⁱ⁺¹and ⁱ⁺¹T_ζ ⁱare right and left pseudo inverses of each other, respectively.
As done above, it is assumed that all branches, active as well as inactive, are assigned to their respective networks and that this particular branch order is fixed for both topologies. Also, for convenience of derivation, the state variables are selected to be inductor currents and capacitor voltages; whereas the dimensions of respective state vectors are kept constant by utilizing (5.36)-(5.42). Then the state variables for the instance of commutation between two topologies can be related as
i _x ⁱ⁺¹ =B _i+1 ^base i _i ^L =B _i+1 ^base(B _i ^L)^T i _x ⁱ=ⁱ⁺¹ T _L ⁱ i _x ⁱ (5.51)
v _y ⁱ⁺¹ =A _i+1 ^base v _i ^C =A _i+1 ^base(A _i ^C)^T v _y ⁱ=ⁱ⁺¹ T _C ⁱ v _y ⁱ (5.52)
If the appropriate redundant states are included in i_x ⁱ, i_x ⁱ⁺¹and v_y ⁱ, v_y ⁱ⁺¹, then the transformation matrices ⁱ⁺¹T_L ⁱand ⁱ⁺¹T_C ⁱ, under conditions of any proper commutation, can always be made to have full rank. Thereafter, ⁱT_ζ ⁱ⁺¹and ⁱ⁺¹T_ζ ⁱbecome inverses of each other.
The second step in deriving an algorithm for global state space continuity consists of applying results of the first step recursively to all topologies of the sequence (5.1). During the time interval of topology s_i, network N=(G, P, s_i) is modeled using 4 transformation of variables of the form
{tilde over (x)}_ζ=⁰K_ζ ⁱx_ζ ⁱ (5.53)
where for ζ=L, C; x_ζ ⁱdenotes i_x ⁱor v_y ⁱ, and x _ζ is the respective state vector transformed by ⁰K_ζ ⁱsuch that state variables are continuous starting from initial topology so all the way up to topology s_i. This step represents an inductive assumption. Then, in order to keep continuity of {tilde over (x)}_ζ during the switching from s_ito s_i+1, based on (5.49), the following relation must be satisfied
{tilde over (x)} _ζ=⁰ K _ζ ⁱⁱ T _ζ ⁱ⁺¹ x _ζ ⁱ⁺¹ (5.54)
From (5.54), it follows that the transformation of variables for the topology s_i+1is in fact
$\begin{matrix} {\overline{x}}_{ζ} = {}^{0}K_{ζ}^{i + 1} x_{ζ}^{i + 1} where & (5.55) \\ {}^{0}K_{ζ}^{i + 1} = \prod_{n = 0}^{n = i} {}^{n}T_{ζ}^{n + 1} & (5.56) \end{matrix}$
Based on (5.50), and very similar to (5.54), it is clear that the inverse of ⁰K_ζ ⁱ⁺¹can also be found as
$\begin{matrix} {}^{i + 1}K_{ζ}^{0} = {}^{i + 1}T_{ζ}^{i} {}^{i}K_{ζ}^{0} = \prod_{n = 0}^{n = i} {}^{i + 1 - n}T_{ζ}^{i - n} & (5.57) \end{matrix}$
The transformation matrices (5.56)-(5.57) can be used in (5.42) in order to achieve {tilde over (x)}_ζ∈C. The same transformation of variables can be applied to the more general system of DAEs (4.10), (4.21), and (4.45)-(4.50), in order to obtain global continuity of state variables.

Implementation

The state selection method discussed herein provides a convenient structure of DAEs for the ASMG. This methodology can be readily implemented on a computer using a matrix representation of a graph corresponding to the global switched network N=(G, P, S). Thereafter, all DAEs for elementary networks can be readily assembled based on the TCF for each topology of N=(G, P, S). However, a brute-force implementation of the ASMG involving sparse topological matrices together with sparse parameter matrices would result in very long simulation run times in the case of electrical networks with time-varying and/or non-linear parameters and a large number of branches. Therefore, an efficient numerical implementation of the ASMG-generated DAEs is a very important issue that will now be addressed.
Instead of employing matrices for the representation of graphs in N=(G, P, S), it is possible to use a collection of linked lists or arrays with multiple indirections. Running spanning tree algorithms on the resulting, more compact data structures requires less CPU time. Utilizing this data structure, the network identification procedure produces arrays containing tree and link branches for each elementary network. Similar arrays can also be assembled by reading the indices of the appropriate rows and columns of the TCF matrix (3.48) or (3.112). Thereafter, it is possible to obtain arrays of branches forming a cycle corresponding to each link branch. Each such array represents a set of branches forming a loop headed by its respective link branch. Also, for each network branch, it is possible to assemble an array that records the loop numbers corresponding to the loops in which this branch participates. Such an array represents a loop participation set for the given branch. Similarly, it is possible to form arrays of cutset branches corresponding to each tree-branch, and arrays of cutset participations for all branches. These compact arrays of loop (loop participation) and cutsets (cutset participation) sets can be formed on a network-by-network basis. As it will be shown, assembling DAEs and computing their terms using arrays of branch sets avoids unnecessary operations and, therefore, significantly reduces the computational complexity.

Implementation of State and Output Equations

In order to reduce the computational effort associated with solving the network DAEs, it is instructive to recall those equations in matrix form. In this section, the output equations (4.13), (4.17), (4.36), and (4.44) are considered first. Then attention is turned to the state equations (3.72) and (3.73), and their more general form (4.34) and (4.43).

Output Equations

Following the sequence in which the material was presented in preceding sections, the inductive network is considered first. For convenience, the voltage equation for N_Lwith time-varying parameters is repeated here.
$\begin{matrix} v_{br}^{L} = R_{br}^{L} i_{br}^{L} + \frac{\partial L}{\partial t} {}^{br}i_{br}^{L} + L_{br} \frac{\partial i_{br}^{L}}{\partial L} + e_{br}^{L} & (6.1) \end{matrix}$
The voltage for the k-th branch of N_Lcan be written directly from (6.1) in terms of summations as
$\begin{matrix} v_{br}^{L} (k) = \sum_{l \in B^{L}} R_{br}^{L} (k, l) i_{br}^{L} (l) + \sum_{l \in B^{L}} \frac{\partial L}{\partial t}^{br} (k, l) i_{br}^{L} (l) + \sum_{l \in B^{L}} L_{br} (k, l) \frac{\partial}{\partial t} i_{br}^{L} (l) + e_{br}^{L} (k) & (6.2) \end{matrix}$
where both indices i and k span all branches of N_L. The number of calculations in (6.2) can be significantly reduced if the summations are performed including only those branches that are coupled through mutual resistances and inductances. Since the coupling among the network branches is known ahead of time, it is possible to define sets that compactly store this information. For instance, it is convenient to define M_k ^Ras the set containing indices of all branches that are coupled by a mutual resistance with the k-th branch of N_L. Also, if branch b_khas a non-zero resistance, it is convenient to let M_k ^Rcontain k as well. On the other hand, if b_khas zero resistance, it cannot be coupled resistively to any other branches and consequently M_k ^Ris empty. Similarly, M_k ^Lis defined to contain indices of inductive branches that are coupled with b_kby a mutual inductance. Again, M_k ^Lcontains k only if branch b_kis inductive and is not coupled with other network branches. If b_khas zero inductance, the corresponding set is M_k ^L=Ø. Furthermore, among branches represented in M_k ^Lthere are some that have time-varying inductance. Thus, it is possible to identify a subset M_k ^Li⊂M_k ^Lthat includes all branches with time-varying inductance. Utilizing such sets, (6.2) can be rewritten as
$\begin{matrix} v_{br}^{L} (k) = \sum_{i \in M_{k}^{R}} R_{br}^{L} (k, l) i_{br}^{L} (l) + \sum_{m \in M_{k}^{Li}} \frac{\partial L}{\partial t}^{br} (k, m) i_{br}^{L} (m) + \sum_{n \in M_{k}^{L}} L_{br} (k, n) \frac{\partial}{\partial t} i_{br}^{L} (n) + e_{br}^{L} (k) & (6.3) \end{matrix}$
Assuming that sets M_k ^R, M_k ^L, and M_k ^Liare assembled for each branch of N_L, the vector of branch voltages v_br ^Lcan be computed using (6.3) with reduced expense compared to (6.2).
The same concept can be applied to the capacitive network. The output equation for the time-varying case of N_Cis
$\begin{matrix} i_{br}^{C} = G_{br}^{C} v_{br}^{C} + \frac{\partial C_{br}}{\partial t} v_{br}^{C} + C_{br} \frac{\partial v_{br}^{C}}{\partial t} - j_{br}^{C} & (6.4) \end{matrix}$
Similar sets corresponding to the branches with mutual conductances M_k ^Gand mutual capacitances M_k ^Cand subsets of branches with time-varying capacitances M_k ^Ct ⊂M_k ^Ccan be assembled for the k-th branch of N_C. The expression for computing the k-th component of i_br ^Cbecomes
$\begin{matrix} i_{br}^{C} (k) = \sum_{l \in M_{k}^{C}} G_{br}^{C} (k, l) i_{br}^{L} (l) + \sum_{m \in M_{k}^{Ct}} \frac{\partial C_{br}}{\partial t} (k, m) v_{br}^{C} (m) + \sum_{n \in M_{k}^{C}} C_{br} (k, n) \frac{\partial}{\partial t} v_{br}^{C} (n) - j_{br}^{C} (k) & (6.5) \end{matrix}$
The sets M_k ^R, M_k ^L, M_k ^Li, and M_k ^G, M_k ^C, M_k ^Ctcan be obtained for the initial topology and then updated for each new topology.
This technique of summing over only relevant terms can be extended to the more general form of output equations (4.36) and (4.44). In the case of inductive networks, among branches represented in M_k ^L, there may be some branches for which the inductance depends on network currents. The indices of these branches can be placed in a set M_k ^φ Similarly, among the capacitive branches represented in M_k ^C, there may be some branches for which the capacitance is a function of applied voltage. These branches are collected in M_k ^φ. Clearly, M_k ^φ⊂M_k ^Land M_k ^φ⊂M_k ^C. After obtaining these additional branch sets, components of (4.36) and (4.44) can be computed as follows
$\begin{matrix} v_{br}^{L} (k) = (\sum_{m \in M_{k}^{R}} R_{br}^{L} (k, m) i_{br}^{L} (m)) + (\sum_{m \in M_{k}^{L}} L_{br} (k, n) \frac{\partial}{\partial t} i_{br}^{L} (n)) + (\sum_{i \in M_{k}^{Lt}} \frac{\partial}{\partial t} L_{br} (k, i) i_{br}^{L} (i)) + (\sum_{j \in M_{k}^{ϕ}} \frac{\partial}{\partial i_{br}^{L} (k)} ϕ_{br} (j) \frac{\partial}{\partial t} i_{br}^{L} (j)) + e_{br}^{L} (k) & (6.6) \\ i_{br}^{C} (k) = (\sum_{m \in M_{k}^{G}} G_{br}^{L} (k, m) v_{br}^{C} (m)) + (\sum_{m \in M_{k}^{C}} C_{br} (k, n) \frac{\partial}{\partial t} v_{br}^{C} (n)) + (\sum_{i \in M_{k}^{Ct}} \frac{\partial}{\partial t} C_{br} (k, i) v_{br}^{C} (i)) + (\sum_{j \in M_{k}^{φ}} \frac{\partial}{\partial v_{br}^{C} (k)} φ_{br} (j) \frac{\partial}{\partial t} v_{br}^{C} (j)) + j_{br}^{C} (k) & (6.7) \end{matrix}$

State Equations

For the purposes of numerical implementation, it is convenient to view the state equation, either for N_Lor for N_C, in the following general form
$\begin{matrix} M (x, t) \frac{\partial x}{\partial t} = A (x, t) + g (u, t) & (6.8) \end{matrix}$
where A(x, t) represents all terms that define the state self-dynamics. The forcing term g(u, t) takes into account all external sources and interconnections with other networks. Since all inputs represented by u have the same units, the function g(u, t) will have the form of a summation. Finally, M(x, t) is the mass matrix.
Thus, N_Lwith time-varying parameters is considered first. The state equation for this case can be written as
$\begin{matrix} L_{x} \frac{\partial i_{x}}{\partial t} = - (R_{x} + \frac{\partial L_{x}}{\partial t}) i_{x} - B_{b}^{L} e_{br}^{L} - C_{b}^{LA} v_{br}^{A} - C^{LC} v_{y}^{C} & (6.9) \end{matrix}$
The forcing vector g(u, t), with respect to (6.9), can be expressed as
g ^L(u,t)=−B _b ^L e _br ^L −C _b ^LA v _br ^A −C ^LC v _y ^C (6.10)
Of course, computing (6.10) by the means of full matrix multiplication is not efficient. Instead, it is possible to utilize the topological information about the matrices in (6.10). In particular, recalling KVL matrix (3.50) and expression (3.62), it can be concluded that the k-th component of the forcing vector (6.10) is the sum of external (forcing) voltages taken around the k-th loop with respect to the inductive network. Then, by examining the nonzero entries of the respective KVL matrices B_b ^L, C_b ^LA, and C^LC, it is possible to assemble special sets of loop branches for each term in (6.10). The set of branch indices corresponding to the k-th inductive link-branch is denoted as L_k ^L, which contains the k-th branch. The set L_k ^Le ⊂L_k ^Lincludes only those N_Lloop-branches that have non-zero external voltage sources. Similarly, L_k ^LCcontains the indices of the loop-branches that happen to be in the algebraic network, and L_k ^LCincludes the loop-branches that are in N_C. Clearly, L_k=L_k ^L∪L_k ^LA∪L_k ^LCcontains all loop-branches (complete basic loop) corresponding to the k-th inductive link-branch. These branch sets can be assembled for each new network topology. Based on these compact loop-sets, the k-th component of the forcing vector can be computed with reduced effort as follows
$\begin{matrix} g^{L} (k) = - \sum_{l \in L_{k}^{Le}} B_{b}^{L} (k, l) e_{br}^{L} (l) - \sum_{m \in L_{k}^{LA}} C_{b}^{LA} (k, m) v_{br}^{A} (m) - \sum_{n \in L_{k}^{LC}} C^{LC} (k, n) v_{y}^{C} (n) & (6.11) \end{matrix}$
The reduced inductance and reduced resistance matrices in (6.9) are defined as a triple product of appropriate matrices. Using to the typical matrix multiplication process, the entries of these matrices are found as
$\begin{matrix} L_{x} (i, j) = \sum_{m = 1}^{η + h} \sum_{n = 1}^{η + h} B_{b}^{L} (i, m) B_{b}^{L} (j, n) L_{br} (m, n) & (6.12) \\ R_{x} (i, j) = \sum_{m = 1}^{η + h} \sum_{n = 1}^{η + h} B_{b}^{L} (i, m) B_{b}^{L} (j, n) R_{br}^{L} (m, n) & (6.13) \end{matrix}$
The terms under summations in (6.12) and (6.13) are non-zero only when the internal indices m and n correspond to the branches that are members of the i-th and j-th respective loops and to the non-zero entry of the respective parameter matrix. Then, using loop sets earlier denoted as L^L, each (i, j)-th element of L_xand R_xcan be computed with reduced computational effort as
$\begin{matrix} L_{x} (i, j) = \sum_{(m \in L_{i}^{L})} \sum_{(n \in L_{j}^{L})} B_{b}^{L} (i, m) B_{b}^{L} (j, n) L_{br} (m, n) & (6.14) \\ R_{x} (i, j) = \sum_{(m \in L_{i}^{L})} \sum_{(n \in L_{j}^{L})} B_{b}^{L} (i, m) B_{b}^{L} (j, n) R_{br}^{L} (m, n) & (6.15) \end{matrix}$
For networks with time-varying inductances and resistances, (6.14) and (6.15) can be updated by running the summation indices m and n only over those branches that have time-varying parameters.
This method of computing the reduced parameter matrices is advantageous over performing the full triple matrix multiplication. An even better performance can be obtained by using loop participation sets. It can be recalled that in the case of N_Lonly the network loop branches that are linked by inductive links make a contribution to (6.12)-(6.13). Thus, for each k-th branch of N_L, it is possible to assemble a set that contains all network loops (loop indices) in which this branch participates as determined by the local KVL matrix B_b ^L. For convenience, the loop numbering can be made the same as for the rows of B_b ^L. It is useful to define LP_kas the set that stores loop numbers in which the k-th branch takes part. This set can be assembled by recording indices of the non-zero entries of the k-th column of B_b ^L. Thereafter, it is possible to express the contribution of mutual inductance between the m-th and n-th branches in the reduced inductance matrix as follows
$\begin{matrix} Δ L_{x} (i, j) = \sum_{(i \in {LP}_{m})} \sum_{(j \in {LP}_{n})} B_{b}^{L} (i, m) B_{b}^{L} (j, n) L_{br} (m, n) & (6.16) \end{matrix}$
When k=m=n, expression (6.16) defines the contribution of the self inductance corresponding to the k-th network branch. In the same way, the contribution of a mutual resistance between the m-th and n-th branches in the reduced resistance matrix can be expressed as
$\begin{matrix} Δ R_{x} (i, j) = \sum_{(i \in {LP}_{m})} \sum_{(j \in {LP}_{n})} B_{b}^{L} (i, m) B_{b}^{L} (j, n) R_{br}^{L} (m, n) & (6.17) \end{matrix}$
Thus, using (6.16)-(6.17), the reduced parameter matrices in the state equation (6.9) can be computed and updated in a very efficient way.
Considering N_Cwith time-varying parameters, the duality with respect to the previous network can be utilized. The corresponding state equation is
$\begin{matrix} C_{y} \frac{\partial v_{y}}{\partial t} = - (G_{y} + \frac{\partial C_{y}}{\partial t}) v_{y} + A_{a}^{C} j_{br}^{C} - D_{a}^{CA} i_{br}^{A} - D^{LC} i_{x}^{L} & (6.18) \end{matrix}$
where the forcing vector can be expressed as
g ^C(u,t)=A _a ^C j _br ^C −D _a ^CA i _br ^A −D ^LC i _x ^L (6.19)
In this case, it is possible to utilize the topological information available in the TCF (3.48). In particular, recalling KCL expression (3.56), it can be concluded that the k-th component of the forcing vector (6.19) is the sum of currents external to the capacitive network. Moreover, the non-zero entries of the KCL matrices A_a ^C, D_a ^CA, and D^LCcorrespond to some cutset branches. The set of cutset branches corresponding to the k-th capacitive tree branch is denoted as C_k ^C, and the subset C_k ^j ⊂C_k ^Cincludes the branches that have nonzero external current sources. The sets C_k ^CAand C_k ^CLinclude the indices of the cutset branches that are in the algebraic and inductive networks, respectively. Clearly, C_k=C_k ^C∪C_k ^CA∪C_k ^clis the complete set of cutset branches corresponding to the k-th capacitive tree-branch. As before, all branch sets are assembled for each new network topology. Thereafter, the k-th component of the forcing vector can be computed with minimized effort as
$\begin{matrix} g^{C} (k) = - \sum_{l \in C_{k}^{Ct}} A_{a}^{C} (k, l) j_{br}^{C} (l) - \sum_{m \in C_{k}^{CA}} D_{a}^{CA} (k, m) i_{br}^{A} (m) - \sum_{n \in C_{k}^{CL}} D^{LC} (k, n) i_{x}^{L} (n) . & (6.20) \end{matrix}$
The reduced parameter matrices can be computed using cutset participation sets. For convenience, the cutset numbering can be made the same as that for the rows of A_a ^C. The set CP_kincludes the cutset indices in which the k-th capacitive branch participates. By analogy, this set can be assembled by recording indices of the non-zero entries of the k-th column of A_a ^C. Thereafter, the contribution of a mutual capacitance and/or mutual conductance between the m-th and n-th branches in their respective reduced matrices can be expressed as
$\begin{matrix} Δ C_{x} (i, j) = \sum_{(i \in {CP}_{m})} \sum_{(j \in {CP}_{n})} A_{a}^{C} (i, m) A_{a}^{C} (j, n) C_{br} (m, n) & (6.21) \\ Δ G_{x} (i, j) = \sum_{(i \in {CP}_{m})} \sum_{(j \in {CP}_{n})} A_{a}^{C} (i, m) A_{a}^{C} (j, n) G_{br} (m, n) & (6.22) \end{matrix}$
Similarly, when k=m=n, (6.21)-(6.22) define contributions of the self-capacitance and self-conductance of the k-th branch of the capacitive network. Using (6.21)-(6.22), the reduced parameter matrices in (6.18) can be computed and updated efficiently.
Building Algorithm for Networks with Variable Parameters
After the appropriate branch sets have been assembled, the implementation of the inductive and capacitive networks becomes very similar. In this section, a pseudo code implementing the output and state equations is described with respect to the inductive network only. The same techniques and conclusions can be directly applied to the capacitive network.
To implement output equation (6.3), the branch sets earlier defined as M_k ^R, M_k ^L, and M_k ^Ltare used. In order to have convenient and fast access to parameters, it is assumed that all parameter matrices use full matrix storage. If standard matrix-vector multiplication is used to implement (6.2), the computational complexity is
η_6.2(n)=Θ(3n+n) (6.23)
However, since the density of matrices in (6.1) is a function the number of mutual parameters, the complexity of (6.3) is somewhat reduced as compared to Θ(3n²+n). In particular, the complexity of the pseudo code in Code Block 1 implementing (6.3) is
η_6.3(n)=Θ[n+n _s ^R +n _s ^L +n _s ^Lt+2(n _m ^R +n _m ^L +n _m ^Lt)] (6.24)
which, in turn, can be anywhere between Θ(4n) and Θ(3n²+n), and would still have the same upper bound on the order of growth of O(n²). In the previous expressions for complexity, n is the number of branches in the given network, n_s ^R,n_s ^L,n_s ^Ltand n_m ^R,n_m ^L,n_m ^Ltare the total numbers of self and mutual resistive, inductive, and time-varying inductive parameters, respectively. The values of those numbers are readily determined from the branch sets M^R, M^L, and M^Lt.


Code Block 1

	for (∀kεB^L) { /* each k-th branch in N_L */
	/* initialize memory slot with e_br ^L */
	VLnw[k] = EsL[k];
	/* j-th and k-th branches have mutual R */
	for (∀jεM_k ^R) {
	VLnw[k] = VLnw[k] + RLnw[k] [j] * ILnw[j];
	}
	/* j-th and k-th branches have mutual L */
	for (∀jεM_k ^L) {
	VLnw[k] = VLnw[k] + Lnw[k] [j] * dILnw[j];
	}
	/* j-th and k-th branches have mutual L(t) */
	for (∀jεM_k ^Lt) {
	VLnw[k] = VLnw[k] + dLnw[k] [j] * ILnw[j];
	}
	}

Implementation of the state equation (6.8) is accomplished in parts. First, the forcing vector is computed according to (6.10). Second, depending on the type of network parameters and choice of state variables, the vector of state derivatives is calculated. A pseudo code block implementing (6.10) is shown in Code Block 2. The code is similar to that shown in Code Block 1. The computational complexity of the pseudo code in Code Block 2 is
$\begin{matrix} η_{6.10} (n) = Θ (n + \sum_{i = 1}^{n} n_{i}^{Le} + \sum_{i = 1}^{n} n_{i}^{LA} + \sum_{i = 1}^{n} n_{i}^{LC}) & (6.25) \end{matrix}$
which has the same upper bound of O(n²). Here, parameter n denotes the number of inductive link branches or the total number of inductive loops, and n_i ^Leis the number of voltage sources in each i-th inductive loop. Then n_i ^LAand n_i ^LCare the numbers of algebraic and capacitive network branches, respectively, that are also part of i-th inductive loop. These quantities are determined by the branch sets L^Le, L^{LA, and L} _LC.


Code Block 2

	for (k=−1...size(B_x ^L)) { /* each k-th inductive link/loop */
	gL[k] = 0.0; /* initialize memory slot */
	/* j-th branch in k-th loop has voltage source */
	for (∀jεL_k ^Le) {
	gL[k] = gL[k] − B[k] [j] * EsL[j];
	}
	/* j-th branch in k-th loop is also in N_A */
	for (∀jεL_k ^LA) {
	gL[k] = gL[k] − CbLA[k] [j] * VAnw[j];
	}
	/* j-th branch in k-th loop is also in N_C */
	for (∀jεL_k ^LC) {
	gL[k] = gL[k] − C_LC[k] [j] * vy[j];
	}
	}

The building algorithm for assembling reduced parameter matrices consists of two parts: initialization and update. The pseudo code implementing these procedures is given in Code Blocks 3-5. Since all reduced parameter matrices can be assembled using their respective loop and loop participation (cutset and cutset participation) branch sets, the building algorithm for the other reduced matrices can be obtained by appropriately modifying the pseudo code given therein. Therefore, only the reduced inductance matrix is discussed.
The computational complexity as well as features of algorithms given in Code Blocks 3-5 are very different. The pseudo code in Code Block 3 can be used to assemble the reduced parameter matrices utilizing loop (cutset) sets as in (6.14). This code may be used to the initialize and update the parameter matrices. The initialization is performed by writing zeros in each memory slot, and then by computing the appropriate contribution due to the constant inductances into an auxiliary storage matrix L_x. Pre-computing the contribution of all time-invariant inductances into a separate matrix suggests this approach
L _x(t)=L _x +ΔL _x(t) (6.26)


Code Block 3

	for (i = 1...size(B_x ^L)) { /* i-th inductive loop */
	for (j = 1...size(B_x ^L)) { /* j-th inductive loop */
	Lx[i] [j] = 0.0; /* initialize memory slot */
	/* m-th branch is in i-th loop */
	for (∀mεL_i ^L) {
	/* n-th branch is in j-th loop */
	for (∀nεL_j ^L) {
	Lx[i] [j] = Lx[i] [j]
	+ BbL[i] [m] * BbL[j] [n] * Lnw[m] [n];
	}
	}
	}
	}


Code Block 4

	/* inductance between m-th and n-th branches of N_L */
	for (∀<m,n>εL_N _L) {
	/* m-th branch is in i-th inductive loop */
	for (∀iεLP_m) {
	/* n-th branch is in j-th inductive loop */
	for (∀jεLP_n) {
	Lx[i] [j] = Lx[i] [j]
	+ BbL[i] [m]BbL[j] [n]L[m] [n];
	}
	}
	}


Code Block 5

	/* variable L between m-th and n-th branches of N_L */
	for (∀<m,n>εL_N _L ^var) {
	/* m-th branch is in i-th inductive loop */
	for (∀iεLP_m) {
	/* n-th branch is in j-th inductive loop */
	for (∀jεLP_n) {
	Lx[i] [j] = Lx[i] [j] + BbL[i] [m] * BbL[j] [n]
	*(Lvar[m] [n] − Lold[m] [n]);
	}
	}
	/* record used value for the future call */
	Lold[m] [n] = Lvar[m] [n];
	}

where ΔL_x(t) is the respective contribution due to remaining variable inductances. Then the procedure for updating L_x(t) can be implemented by initializing each memory slot of L_x(t) with pre-computed values from corresponding slots in L_x. This technique, however, has the disadvantage of having to copy the entire matrix L_xat each update. In large networks where only few parameters are variable, copying the entire matrix L_xover to L_x(t) may represent significant and unnecessary overhead.

Utilizing the loop participation sets according to (6.16)-(6.17), the reduced parameter matrices can be assembled as shown in Code Block 4, without the disadvantage of having to copy the entire matrix for each update. Thus, the code in Code Block 4 avoids the scheme (6.26) altogether. However, if the code in Code Block 4 is used to update the reduced matrix, the original matrix will be destroyed which, in turn, makes it difficult to carry the contribution due to time-invariant parameters from one update to the next. In the algorithm shown in Code Block 3, this function is accomplished by copying the entire matrix; whereas Code Block 4 updates only the relevant entries. One way of using the algorithm in Code Block 4 for systems with both time-varying and time-invariant parameters is to “undo” the previous update before performing a new one. This two-step update procedure is equivalent to updating the matrix once with the difference between the old and the new values of the variable parameters. The pseudo code illustrating the update procedure is given in Code Block 5. Therein, the reduced inductance matrix is updated due only to changes in variable inductances. The implementation of this update algorithm requires auxiliary static memory for storing previous update values. This storage should be of the same size and type as the original reduced parameter matrix. It is also noted that the subtraction of the old and new inductances should be performed first as indicated by the parentheses so as to reduce the round-off errors due to finite precision machine arithmetic.
The complexity of the algorithms in Code Blocks 3-5 is a function of the size of the loop and loop participation sets. The sizes of these sets and their structures are, in turn, determined by the self parameters and the mutual coupling between network branches. Therefore, the complexity of algorithms Code Blocks 3-5 is highly system-dependent. Analyzing each algorithm for the worst case, which corresponds to all branches having variable parameters with all branches coupled with each other, results in similar complexity for all algorithms and reveals little about their actual performances with respect to practical networks. Instead, the problem of assembling the reduced parameter matrices using algorithms in Code Block 3 and Code Block 4 should be considered with respect to some typical cases. Expressions for the complexities can be significantly simplified using certain assumptions. Thus, for the purpose of derivation, it is assumed that there is no mutual coupling. Then, inspecting Code Block 3 for implementing (6.14) with this assumption in mind, the complexity can be expressed as
η_6.14(n)=Θ[n ²( m ²+1)] (6.27)
where n is the number of inductive loops, and m is the average number of branches in each loop. Denoting m_ito be the number of branches in the i-th loop, (6.27) can be developed further as
$\begin{matrix} η_{6.14} (n) = Θ [{n^{2} (\frac{1}{n} \sum_{i = 1}^{n} m_{i})}^{2} + n^{2}] = Θ [{(\sum_{i = 1}^{n} m_{i})}^{2} + n^{2}] & (6.28) \end{matrix}$
With respect to (6.16), the complexity of the code in Code Block 4 is determined by the sizes of the loop participation sets LP. In particular, assuming only self inductances, the computational effort for the contribution from one branch is proportional to the square of the number of loops in which this branch participates. Denoting the number of loops in which the i-th branch participates as p_i, including initialization, the expression for the complexity for the code in Code Block 4 becomes
$\begin{matrix} η_{6.16} (n) = Θ (\sum_{i = 1}^{n} p_{i}^{2} + n^{2}) & (6.29) \end{matrix}$
Both (6.28) and (6.29) demonstrate significant computational savings when compared to Θ[n²(n_br ²+1)], which is the cost of triple matrix multiplication. Here, n_brdenotes the total number of branches in the given network. A closer look at (6.28) and (6.29) suggests that for large networks, the advantage of the building algorithm presented here should become more noticeable.
The run-time computational routines in one embodiment of the present invention will now be discussed in relation to FIG. 5. ODE solver 41 maintains state vector x, which is provided to inductive links current calculator 43 (for state variables relating to inductive elements) and capacitive trees voltage calculator 45 (the portions x_Crepresent state variables relating to capacitive elements). Inductive links current calculator 43 calculates the current i_link ^Lin the inductive link branches of the circuit, providing that information to resistive network algebraic equation calculator 47 and inductive network state/output equation evaluator 49. Likewise, capacitive trees voltage calculator 45 calculates the voltages v_tree ^Cin capacitive tree branches as discussed above, and provides that information to resistive network algebraic equation component 47 and capacitive network state/output equation evaluator 51.
Resistive network algebraic equation evaluator 47 uses i_link ^Land v_tree ^Cwith e_br ^Aand j_br ^Ato calculate i_br ^Aand v_br ^A, which are provided to inductive network state/output equation component 49, capacitive network state/output equation component 51, event variable calculator 53, and the system output. Inductive network state/output equation component 49 uses i_link ^Land v_br ^Aalong with inputs e_br ^Land j_br ^Lto determine i_br ^L, v_br ^L, (provided to event variable calculator 53 and the system output) and dx_L/dt (provided to the ODE solver 41). Capacitive network state/output equation component 51 uses i_br ^Aand v_tree ^Calong with e_br ^Cand j_br ^Cto calculate i_br ^C, v_br ^C, (provided to event variable calculator 53 and the system output) and dx_C/dt (provided to the ODE solver 41).
The branch voltages and currents are output as vectors i_brand v_br. These values are used with u_brby event variable calculator 53 to produce event variable z_xnp, which is passed to the ODE solver 41. In addition to solving the differential equations for the system, ODE solver 41 monitors for negative-to-positive zero crossings of Z_xnp. If a zero crossing is encountered, which indicates that a switch or switches are opening or closing, the state selection algorithm (part of state model generator 31 in FIG. 3) if invoked to establish a new set of state variables and to update the branch sets and matrices used by the state equation building algorithm discussed above.
Inductive links current calculator 43 will now be discussed in relation to FIG. 6. At decision block 61, it is determined whether currents or fluxes have been used to define the state variables for the system. If currents have been used, then the state variables from x_Lare already in the proper form and are provided as output. Otherwise, if fluxes have been used, it is determined at decision block 63 whether constant or variable inductance parameters are present. If the parameters are constant, the output is calculated at block 67 as L_x ⁻¹x_L. Alternatively, if the inductance parameters are variable, block 69 determines the output vector by solving L_xi_link ^L=x_L.
The analogous calculations for capacitive tree voltage calculator 45 will now be discussed in relation to FIG. 7. At decision block 71, it is determined whether voltages or charges have been used to define the state variables for the system. If voltages have been used, then the state variables from x_Care already in the proper form and are provided as output. Otherwise, if charges have been used, it is determined at decision block 75 whether constant or variable capacitance parameters are present. If the parameters are constant, the output is calculated at block 77 as C_y ⁻¹x_C. Alternatively, if the capacitance parameters are variable, block 79 determines the output vector by solving C_yv_tree ^C=x_C.
Resistive network algebraic calculation block 47 will now be discussed in more detail in relation to FIG. 8. At decision block 81, it is determined whether the resistive network parameters are constant or variable. If they are constant, the output current vector is calculated as shown in block 83. If it is determined at decision block 81 that the resistive network parameters are variable, the current vector for the algebraic network branches is determined in block 85 by solving the equation shown therein. After the current vector is determined in either block 83 or block 85, the voltage vector is determined at block 87, as shown therein.
The inductive network state and algebraic equation determining block 49 is shown in more detail and will now be discussed in relation to FIG. 9. First, at block 90, the inductive network branch currents i_br ^Land forcing term g^Lare computed. At decision block 91 it is determined whether the inductive network contains constant or variable parameters. If the parameters are constant, then the type of state variables being used is identified at decision block 92. For currents, the derivatives with respect to time of the inductive portion of the state vector and the currents in the inductive link branches are calculated. Similarly, for state variables defined as fluxes, those values are calculated at block 94. In either event (from either block 93 or block 94), the time-derivative of inductive branch currents and the voltages for inductive branches are determined at block 95, and the outputs are generated.
If it is determined at decision block 91 that variable parameters are present, then it is determined at decision block 96 whether currents or fluxes have been selected as state variables for the inductive portion of the network. If currents have been selected, then the time-derivatives of X_Land the inductive link-branch currents are determined at block 97. Correspondingly, if fluxes are used as state variables, then those values are determined at block 98. In either event (from either block 97 or block 98), the time-derivative of inductive branch currents, as well as the inductive branch voltages, are determined at block 99, and the output of overall block 49 is generated.
Analogously to other discussion in relation to FIG. 9, the capacitive network state in algebraic equation component 51 will now be discussed in relation to FIG. 10. At block 100, values are calculated for the voltages v_br ^Cand forcing term g^Cfor the capacitive network. At decision block 101, it is determined whether constant or variable parameters are present. If constant parameters are being used, then it is determined at decision block 102, whether voltages or charges are used as state variables for the capacitive network. In the former case, the time-derivative of the state variables relating to the capacitive network and the voltages for the capacitive tree branches are calculated at block 103. Likewise, if charges are being used as state variables, those values are calculated at block 104. In either event (whether from block 103 or block 104), the system calculates the capacitive branch currents and the time-derivative of the capacitive branch voltages at block 105, and the output values are generated.
If at decision block 101 it is determined the variable parameters are present, then the type of state variables (voltages or charges) is determined at decision block 106. If voltages are being used, then the time-derivatives of the state variables relating to capacitive branches, as well as the voltages relating to capacitive tree branches, are determined at block 107. If charges are being used for state variables, then the equations of block 108 are solved. In either event (whether from block 107 or block 108), the capacitive branch currents and the time-derivative of the capacitive branch voltages are determined at block 109, and the output of overall block 51 is generated.

Built-in Switching Logic

Although, as shown in FIG. 4, the ASMG has a single switch branch, different logic may be specified by the user to determine when a given switch is opened or closed. Four built-in switch types, each implementing specific switching logic, have been considered in an exemplary embodiment of the invention. These switch types were selected to represent many common solid-state switching devices such as diodes, thyristors, transistors (MOSFET, BJT, IGBT, for example), triacs, and the like. In general, the built-in switching logic does not permit the opening of switches that would cause discontinuities of currents in inductors and/or current sources, as well as closing of switches that would cause discontinuities of capacitor voltages and/or voltage sources. Such attempts would violate Kirchhoff's current law and/or Kirchhoff's voltage law and, therefore, are not allowed. In some embodiments, violation of KCL, KVL, and/or energy conservation principles results in appropriate error messages. The four types of switches built into the system will now be discussed in relation to FIG. 11, in which the switching control signal is denoted as u_br.

Unlatched Bidirectional Switch (UBS)

The switch of this type can conduct current in either direction when the switching control signal u_bris greater than zero, and block positive or negative voltages otherwise. That is, u_br>0 implies v_br=0, while u_br≦0 implies i_br=0. The switch can be opened or closed at any instant of time by controlling variable u_br, subject to KCL, KVL, and energy conservation principles.

Unlatched Unidirectional Switch (UUS)

The switch of this type is similar to the UBS, but can conduct current only in the positive direction. That is, u_br>0 and i_br>0 implies u_br=0, while u_br≦0 implies i_br=0. In any case, i_br≧0. The switch is closed when the variable e₁=min(u_br,v_br) crosses zero going negative-to-positive. The switch becomes open when the variable e₂=min(u_br,i_br) crosses zero going positive-to-negative. This switch should not be used to short-circuit a loop of capacitors and/or voltage sources, unless the resulted sum of the voltages equals to zero.

Latched Bidirectional Switch (LBS)

The switch of this type is also similar to the UBS with the exception that it can be open only at the current-zero-crossing. That is, the switch is closed when u_br>0, which results in u_br=0. In this case, the event variable e₁=u_bris monitored for negative-to-positive zero crossing. In order to open, the switch state controlling variable u_brhas to become negative first. After that, the switch automatically opens at the next current zero crossing. This logic can be used to represent AC arcing switch or an ideal solid-state triac.

Latched Unidirectional Switch (LUS)

The switch of this type is similar to the LBS with the exception that it can conduct current only in positive direction. That is, u_br>0 and i_br>0 implies u_br=0. The switch is closed when the variable e₁=min(u_br, v_br) crosses zero going negative-to-positive. The switch becomes open when the variable e₂=min(u_br,i_br) crosses zero going positive-to-negative. This switch can be used to model an ideal thyristor. These four switch types can be advantageously integrated into the circuit simulation system routines. For example, switching analysis and topology evaluation for state selection can be optimized using the additional information inherent in each switch type, as will occur to those skilled in the art.
Many variations on the above-disclosed exemplary embodiments will be apparent to those skilled in the art. For example, the processes discussed above can be implemented on many different hardware platforms running many different operating systems. The computational portions might be implemented as an integrated portion of the overall system, or might rely on a computational platform such as ACSL (published by The Aegis Technologies Group, Inc., of Huntsville, Ala., USA), or MATLAB/Simulink (published by The MathWorks, Inc., of Natick, Mass., USA).
Furthermore, in various embodiments, the processes might be carried out on a single processor or be distributed among multiple processors.
Still further, the time-domain steps in the simulation might be the same throughout the system, or might vary between variables, for a single variable (over time), between portions of the circuit being simulated, or other divisions as would occur to one skilled in the art. Likewise, the numerical techniques used to perform integration and/or differentiation (e.g., trapezoidal, NDF, or Adams techniques), the rates and maximum error parameters can also vary and might be consistent across the system, vary among portions of the circuit, change over time, or otherwise be applied as would occur to one skilled in the art. Various embodiments of the present invention will also use different techniques to revise the state equations stored for each topology. In some cases, the data structure(s) that describe the state equations before a topology change event are modified only as much as necessary to take into account the new topology (i.e., only the changed portions of the circuit). Alternatively or additionally, new state equations may be derived in whole or in part for one or more topologies. In still other embodiments, a cache is maintained of state equations for some or all of the topologies that are encountered. A wide variety of caching strategies and techniques are available for use with the present invention, as would occur to one skilled in the art.
The parameters of the branches in the system can also be updated during the simulation, using a variety of strategies. In some embodiments, a data structure reflecting the constant parameters is maintained. Each time the variable parameter values are updated in the system, the constants are copied into a new data structure, and the variable parameters are added to them. In other embodiments, the parameters for the circuit at a given time t_iare stored in a data structure. When the parameters are being updated for processing time t_i+1, the variable parameters from time t_iare subtracted from the values in the data structure, then updated for time t_i+1, and the new values are added to the values in the data structure.
While a detailed description of one form of nodal analysis is presented herein, many additional and alternative techniques, approaches, and concepts may be applied in the context of the present invention. Likewise, a variety of different circuits can be simulated using the present system, such as all-inductive circuits with no capacitors, circuits with resistors and voltage or current sources, but neither inductors nor capacitors, to name just a few examples. Parameters may be set for a branch or combination of branches to model ideal or real circuit components or sub-circuits, such as diodes, transistors, thyristors, or triacs.
Using the techniques discussed above for time-varying parameters, multiple models can be used for a single physical component or sub-circuit. For example, a detailed, computationally intensive model might be used for a component when it has a rapidly varying input. Then, when the input has settled to a slower-varying state, a simpler, less computationally intensive model may be substituted for the complex one.
Yet further, the data structures used to represent information in various embodiments of the present invention vary widely as well. The stated structures may be optimized for programming simplicity, code size, storage demands, computational efficiency, cross-platform transferability, or other considerations as would occur to one skilled in the art.
Partitioning of branch sets as discussed herein may employ a wide variety of algorithms as would occur to one skilled in the art. Some spanning tree algorithms that are well adapted for use with the present invention are presented in T. H. Cormen, C. E. Leiserson, R. L. Rivest, Introduction to Algorithms, MIT Press, McGraw Hill, 1993; and R. E. Tarjan, Data Structures and Network Algorithms, Bell Laboratories, Murray Hill, 1983. The implementation of partitioning algorithms to the data structures involved will likely be a consideration for each particular implementation of the present invention.
All publications, prior applications, and other documents cited herein are hereby incorporated by reference in their entirety as if each had been individually incorporated by reference and fully set forth.
While the invention has been illustrated and described in detail in the drawings and foregoing description, the same is to be considered as illustrative and not restrictive in character, it being understood that only the preferred embodiments have been shown and described and that all changes and modifications that would occur to one skilled in the relevant art are desired to be protected. It should also be understood that while the use of the word “preferable,” “preferably,” or “preferred” in the description above indicates that the feature so described may be more desirable in some embodiments, it nonetheless may not be necessary, and embodiments lacking the same may be contemplated as within the scope of the invention, that scope being defined only by the claims that follow. In reading the claims it is intended that when words such as “a,” “an,” “at least one,” “at least a portion,” and the like are used, there is no intention to limit the claim to exactly one such item unless specifically stated to the contrary in the claim.

Claims

1. A method, comprising:

creating one or more data structures sufficient to model an electronic circuit as a collection of n elements consisting of:

zero or more LRV elements, each having at least one of (a) a non-zero inductance parameter L_br, (b) a non-zero resistance parameter r_br, or (c) a non-zero voltage source parameter e_br, but neither a non-zero capacitance parameter, nor a non-zero current source parameter, nor a switch parameter;

zero or more CRI elements, each having at least one of (a) a non-zero capacitance parameter C_br, (b) a non-zero resistance parameter r_br, or (c) a non-zero current source parameter j_br, but neither a non-zero inductance parameter, nor a non-zero voltage source parameter, nor a switch parameter; and

zero or more switching elements, each having a switch state and neither a non-zero inductance parameter, a non-zero capacitance parameter, a non-zero resistance parameter, a non-zero voltage source parameter, nor a non-zero current source parameter; and

automatically generating a first set of state equations from said one or more data structures; and

simulating operation of the electronic circuit by application of said first set of state equations;

wherein n is at least two, and the collection comprises either

an LRV element for which at least two of L_br, r_br, or e_brare non-zero, or

a CRI element for which at least two of C_br, r_br, or j_brare non-zero.

2. The method of claim 1, wherein said simulating comprises producing state output data, the method further comprising:

modifying the parameters in said first set of state equations as a function of said state output data.

3. The method of claim 1, further comprising:

modifying the parameters in said first set of state equations based on a time-varying parameter of at least one element in said collection.

4. The method of claim 1, further comprising:

generating a second set of state equations from said one or more data structures upon the occurrence of a first topology change event.

5. The method of claim 4, wherein said generating said second set of state equations comprises modifying only the subset of said first set of state equations that depend on the one or more switching elements that have changed.

6. The method of claim 4, wherein each unique vector of switch states represents a topology of the overall circuit, and further comprising:

storing said first set of state equations in a cache;

after a second topology change event, determining whether a set of state equations in the cache represents the new topology;

if said determining is answered in the affirmative, using the set of state equations that represents the new topology to simulate operation of the circuit after the second topology change event; and

if said determining is answered in the negative, building a third set of state equations that represents the new topology, and using the third set of state equations to simulate operation of the circuit after the second topology change event.

7. The method of claim 6, further comprising:

storing said second set of state equations in a cache;

after a third topology change event, deciding whether a set of state equations in the cache represents the new topology;

if said deciding is concluded in the affirmative, using the set of state equations from the cache that represents the new topology to simulate operation of the circuit after the third topology change event; and

if said deciding is concluded in the negative, building a new set of state equations that represents the new topology, and using the new set of state equations to simulate operation of the circuit after the third topology change event.

8. A method, comprising:

creating one or more data structures that together store characteristics of a plurality of active branches B^activethat make up a graph of nodes and branches that form a circuit, wherein B^activeconsists of

a set B^Lof zero or more inductive branches, each having a non-zero inductive component but neither a capacitive component nor a variable switch state;

a set B^Cof zero or more capacitive branches, each having a non-zero capacitive component but neither an inductive component nor a variable switch state; and

a set B^Aof additional branches, each having neither an inductive component, nor a capacitive component;

partitioning B^activeinto a first branch set B_tree ^activeand a second branch set B_link ^active, where the branches in B_tree ^activeform a spanning tree over B^active, giving priority in said partitioning to branches not in B^Lover branches in B^L;

sub-partitioning B_link ^activeinto a third branch set B_link ^Land a fourth branch set B_link ^CA, where B_link ^L=B_link ^active∩B^L;

identifying a fifth branch set B^CAas the union of

B_link ^CA,

B^C∩B_tree ^active, and

those branches in B_tree ^activethat form a closed graph when combined with B_link ^CA;

partitioning B^CAinto a sixth branch set {tilde over (B)}_tree ^CAand a seventh branch set {tilde over (B)}_link ^CAwhere the branches in {tilde over (B)}_tree ^CAform a spanning tree over B^CA, giving priority in said partitioning to branches in B^Cover branches not in B^C;

identifying an eighth branch set B_tree ^C={tilde over (B)}_tree ^CA∩B^C;

selecting a set of state variables comprising:

for each branch of B_link ^L, either the inductor current or inductor flux, and

for each branch of B_tree ^Ceither the capacitor voltage or capacitor charge; and

simulating a plurality of states of the circuit using the set of state variables.

9. The method of claim 8, wherein said partitioning steps each comprise an application of a weighted spanning tree algorithm.

10. The method of claim 9 wherein, for some positive numbers w_Land w_C:

for the partitioning of B^active, a minimum spanning tree algorithm is used with weight function

ω_{L} (b_{j}) = {\begin{matrix} w_{L} if branch b_{j} \in B^{L} \\ 0 otherwise \end{matrix}; and

for the partitioning of B^CA, a maximum spanning tree algorithm is used with weight function

ω_{C} (b_{j}) = {\begin{matrix} w_{C} if branch b_{j} \in B^{C} \\ 0 otherwise . \end{matrix}

11. A system, comprising a processor and a computer-readable medium in communication with said processor, said medium containing programming instructions executable by said processor to:

build state equations for a first topology of an electronic circuit having at least two switching elements, wherein each switching element has a switching state;

solve said state equations at time t_ito provide a state output vector, in which at least two elements control the switching states of the switching elements;

calculate the value of a switching variable as a function of the state output vector, wherein the value reflects whether the switching state of at least one of the switching elements is changing; and

if the value of the switching variable at time t_iindicates that at least one of the switching elements is changing, determine a second topology of the electronic circuit for time t_i ⁺ and obtain state equations for the second topology.

12. The system of claim 11, wherein:

said programming instructions comprise a state equation building module, a solver module for ordinary differential equations, and a switching logic module;

said building is performed by the state equation building module;

said solving and calculating are performed by the solver module; and

said determining is performed by the switching logic module.

13. The system of claim 12, wherein said obtaining is performed by said switching logic module.

14. The system of claim 12, wherein said obtaining is performed by said state equation building module.

15. The system of claim 12, wherein:

at a time t_j, at least two switching elements are each either rising-sensitive or falling-sensitive switches, wherein

rising-sensitive switches change switching state if and only if a controlling element of the state vector has passed from a negative value to a non-negative value; and

falling-sensitive switches change switching state if and only if a controlling element of the state vector has passed from a positive value to a non-positive value; and

the function is the arithmetic maximum of

a maximum of all elements of the state vector that control rising-sensitive switches, and

the negative of the minimum of all controlling elements of the state vector that control falling-sensitive switches.

16. A system for simulating electronic circuits, comprising a processor and a computer-readable medium in communication with said processor, said medium containing programming instructions executable by said processor to read element parameters and node connection information from a data stream comprising at least one switch type specification, the at least one switch type specification being selected from the group consisting of:

a unidirectional, unlatched switch;

a bidirectional, unlatched switch;

a unidirectional, latched switch; and

a bidirectional, latched switch; and

wherein said instructions are further executable by said processor automatically to calculate state equations for the circuit given the states of switches specified by said at least one switch type specification.