Initialization and DC Simulation in VHDL 1076.1

Ernst Christen christen@analogy.com


1. Introduction

The purpose of this paper is to define the following aspects related to the simulation of a VHDL 1076.1 model

The table of contents of this document is available here.


2. VHDL 1076.1 Initialization, Quiescent State and Discontinuity Handling


2.1 Domain Signal

[At 2.2-201]

An explicit signal is a signal other that the implicit signals BREAK or DOMAIN, or other than an implicit signal GUARD, or ... (and other modifications in the same paragraph)

[At 4.3-102]

- The implicit signal DOMAIN.

[At 12.6-439]

Finally, the kernel process contains a driver for, and a variable representing the current value of, the implicit DOMAIN signal of type DOMAIN_TYPE.

[There are probably a few more places that must be updated to accommodate DOMAIN.]


2.2 Equations Solved by the Analog Solver

[At 1076.1 LRM 12.4.6-1154]

These characteristic expressions, together with the characteristic expressions from simultaneous statements, form the basic set of characteristic expressions. The set of characteristic expressions processed by the analog solver consists of the basic set and one of five possibly modified augmentation sets. Each augmentation set contains the characteristic expressions corresponding to source quantities and implicit quantities of the form Q'DOT, Q'INTEG, and Q'DELAYED(T) that appear in the text of the model.

[At 1076.1 LRM 12.4.6-1176]

12.4.6.1 Quiescent state augmentation set

Each scalar subelement of each source quantity is a characteristic expression of the quiescent state augmentation set.

Each scalar subelement of each quantity of the form Q'DOT is a characteristic expression of the quiescent state augmentation set.

Each scalar subelement of the prefix Q of each quantity of the form Q'INTEG is a characteristic expression of the quiescent state augmentation set.

The difference between each scalar subelement of each quantity of the form Q'DELAYED(T) and the corresponding scalar subelement of its prefix Q is a characteristic expression of the quiescent state augmentation set.

12.4.6.2 Time domain augmentation set

Each scalar subelement of each source quantity is a characteristic expression of the time domain augmentation set.

The difference between each scalar subelement of each quantity of the form Q'DOT and the time derivative of the corresponding scalar subelement of its prefix Q is a characteristic expression of the time domain augmentation set.

The difference between each scalar subelement of each quantity of the form Q'INTEG and the integral over time from 0.0 to Tc of the corresponding scalar subelement of its prefix Q is a characteristic expression of the time domain augmentation set.

For quantities of the form Q'DELAYED(T) the following rules apply:

12.4.6.3 Discontinuity augmentation set

Each scalar subelement of each source quantity is a characteristic expression of the discontinuity augmentation set.

The difference between each scalar subelement of the prefix Q of each quantity of the form Q'DOT and its value is a characteristic expression of the discontinuity augmentation set.

The difference between each scalar subelement of each quantity of the form Q'INTEG and its value is a characteristic expression of the discontinuity augmentation set.

The difference between each scalar subelement of each quantity of the form Q'DELAYED(T) and its value is a characteristic expression of the discontinuity augmentation set.

12.4.6.4 Frequency domain augmentation set

The difference between each scalar subelement of each spectral source quantity and the value of the expression (ieee.math_complex.exp ( ieee.math_complex.polar_to_complex ( ieee.math_complex.complex_polar' ( magnitude, phase)))), where magnitude and phase denote the corresponding scalar subelements of the magnitude simple expression and phase simple expression from the declaration of the spectral source quantity, is a characteristic expression of the frequency domain augmentation set.

Each scalar subelement of each noise source quantity is a characteristic expression of the frequency domain augmentation set.

The difference between each scalar subelement of each quantity of the form Q'DOT and the corresponding scalar subelement of its prefix Q multiplied by the value of the expression (ieee.math_complex.math_cbase_j * ieee.math_real.math_two_pi * FREQUENCY) is a characteristic expression of the frequency domain augmentation set.

The difference between each scalar subelement of each quantity of the form Q'INTEG and the corresponding scalar subelement of its prefix Q divided by the value of the expression (ieee.math_complex.math_cbase_j * ieee.math_real.math_two_pi * FREQUENCY) is a characteristic expression of the frequency domain augmentation set.

The difference between each scalar subelement of each quantity of the form Q'DELAYED(T) and the corresponding scalar subelement of its prefix Q multiplied by the value of the expression (ieee.math_complex.exp(-T * ieee.math_complex.math_cbase_j * ieee.math_real.math_two_pi * FREQUENCY)) is a characteristic expression of the frequency domain augmentation set.

12.4.6.5 Noise augmentation set

Each scalar subelement of each source quantity is a characteristic expression of the noise augmentation set.

The difference between each scalar subelement of each quantity of the form Q'DOT and the corresponding scalar subelement of its prefix Q multiplied by the value of the expression (ieee.math_complex.math_cbase_j * ieee.math_real.math_two_pi * FREQUENCY) is a characteristic expression of the noise augmentation set.

The difference between each scalar subelement of each quantity of the form Q'INTEG and the corresponding scalar subelement of its prefix Q divided by the value of the expression (ieee.math_complex.math_cbase_j * ieee.math_real.math_two_pi * FREQUENCY) is a characteristic expression of the noise augmentation set.

The difference between each scalar subelement of each quantity of the form Q'DELAYED(T) and the corresponding scalar subelement of its prefix Q multiplied by the value of the expression (ieee.math_complex.exp(-T * ieee.math_complex.math_cbase_j * ieee.math_real.math_two_pi * FREQUENCY)) is a characteristic expression of the noise augmentation set.


2.3 Initialization

[Next is the full initialization, including steps that did not change. Steps are numbered here to allow to reference them in the rationale. They are unnumbered in the LRM.]

[At LRM 12.6.4-634]

  1. existing
    The driving value and the effective value of each explicitly declared signal are computed, and the current value of the signal is set to the effective value. This value is assumed to have been the value of the signal for an infinite length of time prior to the start of simulation.
  2. new
    The value of each explicitly declared quantity is set to its initial value.
  3. modified
    The value of the implicit DOMAIN signal is set to INITIALIZATION_DOMAIN. The value of each implicit signal of the form S'Stable(T) or S'Quiet(T) is set to True. The value of each implicit signal of the form S'Delayed(T) is set to the initial value of its prefix, S. The value of each implicit signal of the form Q'Above(E) is set to the value of the boolean expression Q > E. The value of each implicit quantity of the form T'Reference, T'Contribution, Q'Dot, Q'Integ, Q'Zoh(T, Initial_delay), Q'Ltf(Num, Den) or Q'Ztf(Num, Den, T, Initial_delay) is set to 0.0. The value of each implicit quantity of the form Q'Delayed(T) or Q'Slew(max_rising_slope, max_falling_slope) is set to the initial value of its prefix, Q. The value of each implicit quantity of the form S'Ramp(tr, tf) or S'Slew(max_rising_slope, max_falling_slope) is set to the initial value of its prefix, S.
  4. existing
    The value of each implicit GUARD signal is set to the result of evaluating the corresponding guard expression.
  5. new
    The driver of the implicit break signal is assigned the waveform TRUE after 0 ns.
  6. existing
    Each nonpostponed process in the model is executed until it suspends.
  7. existing
    Each postponed process in the model is executed until it suspends.
  8. new
    The basic set of characteristic expressions is augmented by the quiescent state augmentation set.
  9. new
    The break set is applied to the quiescent state augmentation set.
  10. existing
    The time of the next simulation cycle (which in this case is the first simulation cycle), Tn, is calculated according to the rules of step h of the simulation cycle, below.

Note: The calculation of step 10 will always yield a value of type Universal_Time corresponding to 0 ns.


2.4 Mixed-Mode Simulation Cycle

[The quiescent state is defined as the point where there are no pending events at 0 ns.]

[At 12.6.4-646]

  1. The analog solver is executed. (Tn may be reset as a result.)
  2. The current simulation time, Tc is set to Tn. Simulation is complete when Tn corresponds to TIME'HIGH and there are no active drivers or process resumptions at Tn.
  3. Each active explicit signal in the model is updated. (Events may occur on signals as a result.)
  4. Each implicit signal in the model is updated. (Events may occur on signals as a result.)
  5. If the value returned by the function DOMAIN'EVENT is False, the remainder of this step is skipped. Otherwise, the time domain augmentation set is elaborated and replaces the quiescent state augmentation set.
  6. For each process P, if P is currently sensitive to a signal S and if an event has occurred on S in this simulation cycle, then P resumes.
  7. Each nonpostponed process that has resumed in the current simulation cycle is executed until it suspends.
  8. The time of the next simulation cycle, Tn, is determined by setting it to the earliest of
    1. The value of type Universal_Time corresponding to TIME'HIGH,
    2. The next time at which a driver becomes active, or
    3. The next time at which a process resumes.

    If Tn = Tc, then the next simulation cycle (if any) will be a delta cycle.

  9. If the value of the DOMAIN signal is TIME_DOMAIN, or if Tn = 0.0, the remainder of this step is skipped. Otherwise, the driver of the implicit DOMAIN signal is assigned one of the following waveforms:
    • TIME_DOMAIN after 0 ns if a time domain simulation will follow, or
    • FREQUENCY_DOMAIN after 0 ns if a frequency domain or noise simulation will follow.

    Then, Tn is reset to 0.0. The next simulation cycle will be a delta cycle.

  10. If the next simulation cycle will be a delta cycle, the remainder of this step is skipped. Otherwise, each postponed process that has resumed but has not been executed since its last resumption is executed until it suspends. Then Tn is recalculated according to the rules of step h. It is an error if the execution of any postponed process causes a delta cycle to occur immediately after the current simulation cycle.


2.5 Break Statement

[At 1076.1 LRM 8.14-992]

break_element ::= [ selector_clause ] quantity_name => expression

selector_clause ::= for quantity_name use

A break element without a selector clause is equivalent to a break element with a selector clause whose name denotes the quantity specified in the break element.

In each break element, the name must denote a quantity whose base type is the same as that of the expression and of the named quantity in the selector clause. The quantity name in the selector clause must denote a quantity that is either the prefix Q of a quantity of the form Q'Dot, or that has the form Q'Integ. This quantity is said to be the selector quantity for the break element.

For the execution of a break statement the condition, if present, is first evaluated. A break is indicated if the value of the condition is TRUE or if there is no condition. If a break is indicated, the driver of the implicit BREAK signal is assigned the waveform TRUE after 0 ns and then each break element is evaluated in the order in which the elements appear.

For the evaluation of a break element, the selector clause, the quantity name and the expression are first evaluated. A model is erroneous if it depends on the order of this evaluation. Then the selector quantity, the quantity named in the break element and the value of the expression are added to the break set for the current simulation cycle (see 12.6.4 and 12.6.5.1).


2.6 Handling of Discontinuities

[At 1076.1 LRM 12.6.5-1302]

When the analog solver resumes at time Tc, if the driver of the implicit break signal is active, and if the value of the predefined DOMAIN signal is TIME_DOMAIN, the discontinuity augmentation set is elaborated and replaces the currently active augmentation set. Then the break set is applied to the resulting augmentation set. Next, the analog solver determines an analog solution point based on the assumption that a discontinuity has occurred. Finally the augmentation set that was active when the analog solver resumed at Tc is re-established.

[At 1076.1 LRM 12.6.5-1317]

12.6.5.1 Application of the break set

The application of the break set replaces certain characteristic expressions in the augmentation set.

The difference between each scalar subelement of the named quantity in each break element and the corresponding scalar subelement of the value of the expression is a characteristic expression that replaces one of the following characteristic expressions in the augmentation set.

A break element that has been applied to the augmentation set is removed from the break set.


3. Rationale


3.1 DOMAIN signal

The value of the predefined signal DOMAIN indicates the domain of the analysis:

The DOMAIN signal serves two purposes. It allows a model writer to write behavior that differs in different domains, and it indicates when the quiescent state of the model has been found. DOMAIN is an implicit signal whose driver is owned by the kernel process. As a consequence, it is an error if a process assigns a value to the DOMAIN signal.


3.2 Equations Solved by the Analog Solver

The analog solver's task is to find values for the quantities such that the evaluation of the characteristic expressions yields values close to zero. The characteristic expressions of the model belong to one of two sets: the basic set or an augmentation set.

The basic set contains the characteristic expressions that are a consequence of simultaneous statements, free quantity declarations, branch quantity declarations, terminal associations, and quantity associations. Because the definition of the implicit quantities of the form Q'ZOH, Q'LTF, Q'ZTF, Q'SLEW, S'RAMP, and S'SLEW is based on equivalent blocks their characteristic expressions are also included in the basic set. Characteristic expressions in the basic set are replaced only as a consequence of simultaneous if statements and simultaneous case statements.

Characteristic expressions that are a consequence of the declaration of source quantities or the appearance of quantities of the form Q'DOT, Q'INTEG or Q'DELAYED(T) in the text of a model belong to each of five augmentation sets. Each augmentation set has a very specific area of application:


3.3 Initialization

The purpose of initialization is to give initial values to all dynamic objects in the model, to determine an initial resumption point for each process, and to set up the equations to be solved by the analog solver. The definition modifies some steps of the VHDL 1076 initialization and adds some new steps. It is backward compatible, i.e. if there is no quantity in a model the initialization sequence reduces to the VHDL 1076 initialization sequence, with the exception of step 5, defining the driver of the break signal, which has no backward compatibility consequences because the break signal cannot be accessed in a model.

The following considerations guided the definition of the VHDL 1076.1 initialization sequence:


3.4 Finding the Quiescent State

The VHDL 1076 simulation cycle has been augmented with three additional steps, two of which are specific to finding the quiescent state of the model. The quiescent state is defined as the point where for the first time there are no pending events at 0 ns. This is expressed by step i: if Tn is > 0.0 while DOMAIN = INITIALIZATION_DOMAIN, the driver of the DOMAIN signal is assigned a value that depends of what kind of simulation follows, a time domain or a frequency domain simulation. In a purely digital model this causes an extra delta cycle to occur, otherwise the augmented simulation cycle is backward compatible with the VHDL 1076 simulation cycle.

The analog solver executes at the beginning of the simulation cycle. During a time domain simulation it establishes a sequence of analog solution points between Tc and Tn.. While finding the quiescent state of the model it finds one analog solution point at the Universal_time corresponding to 0.0 during each delta cycle.

Once the DOMAIN signal has been updated the set of equations to be solved by the analog solver is changed: the quiescent state augmentation set is replaced with the time domain augmentation set that constrains implicit quantities in the following way:

The following considerations were made:


3.5 Handling of Discontinuities

When the driver of the break signal is active when the analog solver resumes a discontinuity has occurred. In this case Tn is guaranteed to be equal to Tc. The handling of discontinuities is slightly different when finding the quiescent state compared to during a time domain simulation.

When a discontinuity occurs during a time domain simulation an initial solution at t+ of the discontinuity must be found. Based on the assumption that we are dealing with a physical system the default behavior at a discontinuity is to conserve energy. This is accomplished by forcing, in the discontinuity augmentation set, the values of all quantities of the form Q'INTEG and all quantities Q for which Q'DOT exists in the model to retain their values across the discontinuity. If a discontinuity is to occur on one of these quantities an explicit initial condition for it must be specified with the break statement. The discontinuity augmentation set replaces the time domain augmentation set temporarily, i.e. for finding an analog solution point, when the analog solver resumes and the break signal is active during a time domain simulation.

While finding the quiescent state the quiescent state augmentation set remains in place even at a discontinuity. Energy is not conserved across a discontinuity because the intent is to find the quiescent state quickly. This has some similarity to constraining the time derivatives to 0.0 while finding the quiescent state: both are algorithmic liberties and means to an end.

Initial conditions specified as break elements in a break statement are applied in three different places. First, initial conditions for finding the quiescent point are applied during initialization. Therefore, the analog solver never "sees" the DC equations if such initial conditions have been specified. Then, if a discontinuity occurs while finding the quiescent state the initial conditions after the discontinuity are applied temporarily, while finding one analog solution point. In both these cases the initial conditions specified in the break set replace certain characteristic expressions in the quiescent state augmentation set. The third place of applying initial conditions is at a discontinuity during a time domain simulation. In this case suitable characteristic expressions from the discontinuity augmentation set are replaced by the characteristic expressions formed from the initial conditions.


4. Initialization of DAEs


4.1 Theoretical Background

The solution of a system of DAEs

	F(x, dx/dt, t) = 0				(1)
 

in any continuous interval depends only on the values of the unknowns x and dx/dt at the beginning of the interval [BRENAN]. In particular, this is the case at the beginning of a time domain simulation, i.e. at time t = 0. The values of the unknowns x and dx/dt must satisfy the original system of DAEs at t = t0, and additionally the equations obtained by differentiating (1) w.r.t. time m times:

	dF(x, dx/dt, t0)/dt = 0
 
	d2F(x, dx/dt, t0)/dt2 = 0			(2)
 
	...
 
	d(m)F(x dx/dt, t0)/dt(m) = 0
 

and the set of user defined initial conditions

	B(x, dx/dt, t0) = 0				(3)
 

Here m is the index of the DAE system [BRENAN]. Taking the derivatives of (1) creates new equations, but also new unknowns: higher order derivatives of x. If the user defined initial conditions B are sufficient to determine a unique solution for F, then (1), (2), (3) have a solution (x, dx/dt,... d(m+1)x/dt(m+1)) where x and dx/dt are uniquely determined [BRENAN].

Many of the equations defined by (2) do not actually yield new information and are therefore redundant. A graph-based algorithm is available that determines the minimum set of equations to create [PANTELIDES].

The numerical solution of (1), (2), (3) is very difficult, particularly if the equations are not available in symbolic form. One such algorithm has been published by Leimkuhler et.al. [LEIMKUHLER], but its application is restricted to certain classes of low index DAE systems. No algorithm is known that solves the general initialization problem for an arbitrary system of DAEs.


4.2 DC Operating Point

Circuit simulators traditionally have provided a simpler initialization scheme: DC operating point simulation. The DC operating point is the solution of the equations

	F(x, dx/dt, t0) = 0				(4a)
 
	dx/dt = 0					(4b)
 

The importance of the DC operating point comes from its physical meaning: it represents the quiescent state of the system described by the DAEs, which is a result that a designer is interested in for many systems. Further, the small-signal model of the system, which is used for small-signal analyses (e.g. small-signal AC analysis), is determined by linearizing the large-signal model about the DC operating point. Finally, the DC operating point has traditionally been used in circuit simulators as an initial point for time domain simulations, unless initial conditions were specified by the user.

Using the DC operating point instead of the solution of the general initialization problem as an initial point for time domain simulation in VHDL 1076.1 has both advantages and disadvantages. Important advantages are:

On the other hand, there are some disadvantages.


4.3 User-Specified Initial Conditions


4.3.1 Why Initial Conditions

Users want to specify initial conditions for a variety of reasons:

Some of these reasons are questionable, particularly in the context of a DC operating point, because they are likely to conflict with the dx/dt = 0 property of a DC operating point.


4.3.2 Forms of Initial Conditions

The form initial conditions take depends on the application. In conservative systems it is common to have initial conditions specifying the values of branch quantities. It is rare to have initial conditions that involve several branch quantities. A similar statement can be made for free quantities in signal flow system models. However, in more advanced control system problems initial conditions relating several quantities are common.

The specification of a particular initial condition in a design entity may be conditional depending on whether a generic has been given a value other than the default value.


4.4 Multiple Solutions

The equations to be solved during initialization can have more than one solution if the problem is not convex. An example is a simple flip-flop, which has three solutions: two are stable and one is unstable. The problem of multiple solutions has been studied for DC operating points in electrical systems [GREEN], and methods have been proposed that allow to find all such solutions. Here we are more interested in specifying that a particular solution be found.

There is no theory available on this topic. Circuit simulators have traditionally used a combination of local and global methods to precondition the solution in a particular way. An example of a local method is the OFF parameter on an instance of a SPICE semiconductor device model: its effect to constrain certain currents to zero is local to one instance. An example of a global method is the .nodeset command in SPICE, which roughly allows to specify approximate values for node voltages.


5. Alternatives Considered


5.1 Boolean Signal QUIESCENT and Impure Function DOMAIN

Keeping these things apart made the simulation cycle slightly simpler, because there was no need to indicate what simulation followed quiescent state. However, there was some overlap in functionality and QUIESCENT was True during time domain simulations, which didn't look right.


5.2 QUIESCENT Indicates no A/D and D/A Interaction

This is an alternative to section 3.4. Rather than no pending events at 0, this version quiescent point as no A/D and D/A interaction.

[Before 12.6.4-646: The VHDL simulation cycle is modified to include the following steps prior to step a)]

w) The driver of the implicit QUIESCENT signal is assigned the waveform True after 0 ns if its effective value is False and the driver of the BREAK signal is inactive.

x) The analog solver is executed. (Tn may be reset as a result.)

y) The driver of the implicit QUIESCENT signal is assigned the waveform True after 0 ns if its effective value is False and the drivers of all implicit signals of the form Q'Above(E) are inactive.

[Step f) is modified to correctly handle detection of the quiescent state]

f) 4) If the value of the QUIESCENT signal is False, the value of type Universal_time corresponding to 0 ns.


5.3 Dominique Rodriguez' proposal for initialization

1. The driving value and the effective value of each explicitly declared signal are computed, and the current value of the signal is set to the effective value. This value is assumed to have been the value of the signal for an infinite length of time prior to the start of simulation.

2. The value of each explicitly declared quantity is set to its default value.

3. The value of each implicit quantity of the form Q'DOT and Q'INTEG is set to 0.0; the value of each implicit quantity of the form Q'DELAYED(T) is set to the initial value of Q.

4. The value of each implicit signal of the form S'Stable(T) or S'Quiet(T) is set to True. The value of each implicit signal of the form S'Delayed(T) is set to the initial value of its prefix, S. The value of each implicit signal of the form Q'Above(E) is set to the value of the corresponding boolean expression.

5. The value of each implicit GUARD signal is set to the result of evaluating the corresponding guard expression.

6. The value of the implicit BREAK signal is set to the projected output waveform "TRUE after 0 ns".

7. Each nonpostponed process in the model is executed until it suspends.

8. Each postponed process in the model is executed until it suspends.

9. goto step 16 if no analog statements are present or if the first entry of the projected output waveform of each driver in the analog interface set is identical to the value of the respective drivers on the previous iteration of the cycle (two cycles minimum).

10. all the signals are updated.

11. The analog solver is executed.

12. The implicit signals of the form Q'Above(E) are updated.

13. Discard the event queue (to be defined, this means that all the driver of all the signal are emptied except the current value)

14. RE-Elaborate the process statements, (in order to re-run from the first statement), all the signal declarations (to be resolved <= all the shared variable declarations)

15. initialization continue in step 4.

16. The time of the next simulation cycle (which in this case is the first simulation cycle), T_n, is calculated according to the rules of step f of the simulation cycle, below.


5.3.1 Dominique's Comments

It has been noticed that the part containing the statements that appear between the BEGIN keyword and the first WAIT statement of a process can be called the init part of the process. This is a conceptual view, because the first WAIT statement can be non-unique. However, This part (which can be empty) is commonly used to set the internal state of the device (ie a finite state machine).

In a mixed device, some quantities of the digital interface set can be used to set this internal initial state. This is why it is important to ensure that the init part of the process is executed with the correct value of all the quantities of the digital interface set.

At the first execution of all the processes, the quantities of the digital interface set have their default value, and the first execution of the analog solver can change this values. This is why it is proposed to re-execute the init part of the process with the new value of the quantities of the digital interface set, but also with the original state of the internal variables (OI: the shared variable) and of the signal. It is to notice that the current value of the signals have not changed since no delta has been performed. However, the driver must be reset in order to empty the event-queue. The last point is to remove all the processes waiting for a timeout from the event queue.

All the processes (normal and post-poned) are re-resumed, all the computation is made using the new-value of the quantities and of the signals of the form Q'Above.

All the processes which use only the value of "digital" signals are going to perform the same elaboration precedently. This mean that the result will be the same for their associated drivers.

The steps 5, 6, 7 and 8 are executed at least twice if there is any simultaneous statements.


5.3.2 Compatibility with digital models

This initialization phase is totally compatible with the VHDL'93 initialization phase. Thus, if no analog objects are defined, no quantities of the digital interface set will be modified as a result of the step 9, so the step 10 will be skipped. Then we have exactly the initialization phase as described in the VHDL'93 LRM.


5.3.3 Properties

For a pure digital model the value of all signals at the end of the initialization phase is the default initial valie. For a purely analog model the analog solver is executed twice (this can easily be fixed). For a mixed-signal model the values of signals at the end of initialization are inconsistent: For example, in an A/D converter all output signals have their default values, only the drivers are loaded. As a consequence, if simultaneous statements depend on signals the small-signal model will be computed wrongly.


5.4 Discontinuity Handling

The problem of re-initializing the DAEs after a discontinuity consists of finding values for all unknowns at t+ of the discontinuity, using the values of the unknowns at t- as starting values. It is similar to initializing the DAEs at time 0, but requires solving the general initialization problem with initial conditions that somehow reflect the discontinuity.

In the literature the problem of handling discontinuities has been approached from different directions, but to our knowledge no general solution has been proposed. Most of the publications deal with linear systems only. Others restrict the problem to a particular class of electrical networks and require knowledge of the element type (e.g. capacitor, switch) of a particular component instance. Although they have shortcomings, two algorithms look promising.

Sincovec et.al [SINCOVEC] proposed an algorithm to find the values of the unknowns at t+ for linear systems. The algorithm has also been applied to nonlinear systems [VLACH], but not much information is available about its performance in the general case. The algorithm does not allow the specification of initial conditions at t+.

Opal and Vlach [OPAL] proposed a method of building a system of equations to be solved to find the values of the unknowns at t+. The method has been proposed for switched electrical networks and requires knowledge of the element type of each component instance.

Both algorithms were not considered further because of their shortcomings.


6. References

[BRENAN]
K.E.Brenan, S.L.Campbell, L.R.Petzold: "Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations", North-Holland, 1989
[PANTELIDES]
C.C.Pantelides: "The Consistent Initialization of Differential-Algebraic Systems", SIAM J.Sci.Stat.Comput., vol.9, no.2, March 1988, pp 213-231
[LEIMKUHLER]
B.Leimkuhler, L.R.Petzold, C.W.Gear: "Approximation Methods for the Consistent Initialization of Differential-Algebraic Equations", SIAM J.Numer.Anal., vol.28, no.1, Feb. 1991, pp 205-226
[GREEN]
M.M.Green, A.N.Wilson: "(Almost) Half of Any Circuit's Operating Points are Unstable", IEEE Trans.Circ.Sys.I, vol.41, no.4, April 1994, pp 286-293
[SINCOVEC]
R.F.Sincovec et.al.: "Analysis of Descriptor Systems Using Numerical Algorithms", IEEE Trans.Aut.Cont., vol.AC-26, no.1, Feb. 1981, pp 139-147
[OPAL]
A.Opal, J.Vlach: "Consistent Initial Conditions of Nonlinear Networks with Switches", IEEE Trans.Circ.Sys., vol.38, no.7, July 1991, pp 698-710
[VLACH]
J.Vlach, J.M.Wojciechowski, A.Opal: "Analysis of Nonlinear Networks with Inconsistent Initial Conditions", IEEE Trans.Circ.Sys.I, vol.42, no.4, Apr.1995, pp 195-200


Revision History

0.1 Jan.13, 1995 Ernst Christen
Original version (incomplete)
0.2 Jan.25, 1995 Ernst Christen
Added initial conditions, streamlined section 4.2., added sections 4.3. and 4.4.
0.3 Feb.24, 1995 Ernst Christen
Included specific comments from LAT meeting
0.4 Jan.11, 1996 Ernst Christen
Major reorganization. Included results on general initialization of DAEs, the algorithm outline for VHDL-A initialization, and information about handling discontinuities.
0.5 Apr.24, 1996 Ernst Christen
Replaced VHDL-A by VHDL 1076.1 or equivalent.
0.6 May 21, 1996 Ernst Christen
Included definitions and rationale for initialization phase. Updated section on discontinuities. Changed specification of initial conditions to use BREAK statement.
0.7 July 29, 1996 Ernst Christen
Defined QUIESCENT signal. Updated initialization and simulation cycle based on feedback from last meeting and including QUIESCENT.
0.8 Sept.6, 1996 Ernst Christen
Added alternative simulation cycle where QUIESCENT becomes true immediately before time is about to advance for the first time.
0.9 Nov. 15, 1996 Ernst Christen
Created new section on alternatives. Moved original proposal to alternatives. Included Dominique Rodriguez' proposal in alternatives.
1.0 Dec. 10 1996 Ernst Christen
Merged QUIESCENT signal and DOMAIN function into DOMAIN signal.
1.1 March 10, 1997 Ernst Christen
Reorganized document. Introduced definitions for basic set of characteristic expressions and five augmentation sets. Updated initialization sequence to reflect new implicit quantities and augmentation sets. Updated simulation cycle to include changing equations when the quiescent state has been reached. Provided a definition for the break element in a break statement that includes a mechanism to select which characteristic expression is to be replaced. Formalized the handling of discontinuities. Updated rationale to reflect new developments.