Author: Pascal Denoyelle, September 22, 2005

You may refer to the source code of this example, found here.

There is also a PDF version of this document, viewable online here.

In this sample, a LC oscillator initialized with a given voltage across the capacitor and a null current through the inductor provides the energy to a load resistance through a half-wave rectifier consisting of an ideal diode (see fig 1: Electrical oscillator with half-wave rectifier).

Only the positive wave of the oscillating voltage across the LC is provided to the resistor. The energy is dissipated in the resistor resulting in a damped oscillation.

This type of non-smooth dynamical system consists of :

- a time invariant linear dynamical system (the oscillator). The state variable of this system is denoted by \(x\).
- a non-smooth law describing the behaviour of the diode as a complementarity condition between current and reverse voltage (variables (\(y,\lambda\)) )
- a linear time invariant relation between the state variable \(x\) and the non-smooth law variables (\(y,\lambda\))

Remark:

In a more general setting, the system’s evolution would be described by a DAE :

\[G \cdot x' = A \cdot x + E \cdot u + b + r\]with \(G , A , E\) matrices constant over time (time invariant system), \(u , b\) source terms functions of time and \(r\), a term coming from the non-smooth law variables : \(r = B \cdot \lambda + a\) with \(B , a\) constant over time. We will consider here the case of an ordinary differential equation :

\[x' = A \cdot x + E \cdot u + b + r\]and an initial value problem for which the boundary conditions are \(t_0 \in \mathbb{R} , x(t_0)= x_0\).

In the linear time invariant framework, the non-smooth law acts on the linear dynamical system evolution through the variable \(r = B \cdot \lambda + a\). Reciprocally, the state variable \(x\) acts on the non-smooth law through the relation \(y = C \cdot x + D \cdot \lambda + F \cdot u + e\) with \(C , D , F , e\) constant over time.

It is a complementarity condition between y and \(\lambda\) : \(0 \leq y \, \perp \, \lambda \geq 0\). This corresponds to the behaviour of the rectifying diode, as described in Non Smooth laws.

The equations come from the following physical laws :

- the Kirchhoff current law (KCL) establishes that the sum of the currents arriving at a node is zero,
- the Kirchhoff voltage law (KVL) establishes that the sum of the voltage drops in a loop is zero,
- the branch constitutive equations define the relation between the current through a bipolar device and the voltage across it

Refering to fig 1: Electrical oscillator with half-wave rectifier, the Kirchhoff laws could be written as:

\[\begin{split}\begin{array}{l}
v_L = v_C\\
v_R + v_D = v_C\\
i_C + i_L + i_R = 0\\
i_R = i_D
\end{array}\end{split}\]

while the branch constitutive equations for linear devices are:

\[\begin{split}\begin{array}{l}
i_C = C v_C'\\
v_L = L i_L'\\
v_R = R i_R
\end{array}\end{split}\]

and last the “branch constitutive equation” of the ideal diode that is no more an equation but instead a complementarity condition:

\[0 \leq i_D \, \perp \, -v_D \geq 0\]

This is illustrated in fig 2: Non-smooth and smooth characteristics of a diode, where the left-hand sketch displays the ideal diode characteristic and the right-hand sketch displays the usual exponential characteristic as stated by Shockley’s law.

After rearranging the previous equations, we obtain:

\[\begin{split}\left( \begin{array}{c}
v_L'\\
i_L'
\end{array} \right)
=
\left( \begin{array}{cc}
0 & \frac{-1}{C}\\
\frac{1}{L} & 0
\end{array} \right)
\cdot
\left( \begin{array}{c}
v_L\\
i_L
\end{array} \right)
+
\left( \begin{array}{c}
\frac{-1}{C}\\
0
\end{array} \right)
\cdot i_D\end{split}\]

that fits in the frame of ref{sec-def-NSDS} with

\[\begin{split}x =
\left( \begin{array}{c}
v_L\\
i_L
\end{array} \right)\end{split}\]

and

\[\lambda = i_D\]

We recall that the \(r = B \cdot \lambda + a\) equation is expressed with

\[\begin{split}r =
\left( \begin{array}{c}
\frac{-1}{C}\\
0
\end{array} \right)
\cdot i_D\end{split}\]

from the dynamical equation (Dynamical equation).

Rearranging the initial set of equations yields:

\[\begin{split}-v_D =
\left( \begin{array}{cc}
-1 & 0
\end{array} \right)
\cdot
\left( \begin{array}{c}
v_L\\
i_L
\end{array} \right)
+ R i_D\end{split}\]

as the second equation of the linear time invariant relation with

\[y = -v_D\]

There is just the complementarity condition resulting from the ideal diode characteristic:

\[0 \leq i_D \, \perp \, -v_D \geq 0\]

The integration of the ODE over a time step \([t_i,t_{i+1}]\) of length \(h\) is:

\[\int_{t_i}^{t_{i+1}}x'\,dt = \int_{t_i}^{t_{i+1}}
A \cdot x\,dt + \int_{t_i}^{t_{i+1}}(E \cdot u + b) dt
+ \int_{t_i}^{t_{i+1}}r\,dt\]

The left-hand term is \(x(t_{i+1})-x(t_i)\).

Right-hand terms are approximated this way:

\(\int_{t_i}^{t_{i+1}} A \cdot x\,dt\) is approximated using a \(\theta\)-method

\[\int_{t_i}^{t_{i+1}} A \cdot x\,dt \approx h \theta (A \cdot x(t_{i+1})) + h (1-\theta) (A \cdot x(t_{i}))\]since the second integral comes from independent sources, it can be evaluated with whatever quadrature method, for instance a \(\theta\)-method

\[\int_{t_i}^{t_{i+1}}(E \cdot u + b) dt \approx h \theta (E \cdot u(t_{i+1}) + b(t_{i+1})) + h (1-\theta) (E \cdot u(t_{i}) + b(t_{i}))\]the third integral is approximated like in an implicit Euler integration

\[\int_{t_i}^{t_{i+1}}r\,dt \approx h r(t_{i+1})\]

By replacing the accurate solution \(x(t_i)\) by the approximated value \(x_i\), we get:

\[x_{i+1}-x_i = h \theta (A \cdot x_{i+1}) + h (1-\theta) (A \cdot x_{i}) +
h \theta (E \cdot u(t_{i+1}) + b(t_{i+1})) + h (1-\theta) (E \cdot u(t_{i}) + b(t_{i})) + h r_{i+1}\]

Assuming that \(I - h \theta A\) is invertible, matrix \(W\) is defined as \((I - h \theta A)^{-1}\). We get then:

\[x_{i+1} = W(I + h (1-\theta) A) \cdot x_{i} +
W (h \theta (E \cdot u(t_{i+1}) + b(t_{i+1})) + h (1-\theta) (E \cdot u(t_{i}) + b(t_{i}))) + h W r_{i+1}\]

An intermediate variable \(x_{free}\) related to the smooth part of the system is defined as:

\[x_{free} = W(I + h (1-\theta) A) \cdot x_{i} +
W (h \theta (E \cdot u(t_{i+1}) + b(t_{i+1})) + h (1-\theta) (E \cdot u(t_{i}) + b(t_{i})))\]

Thus the calculus of \(x_{i+1}\) becomes:

\[x_{i+1} = x_{free} + h W r_{i+1}\]

It comes straightforwardly:

\[ \begin{align}\begin{aligned}r_{i+1} =& B \cdot \lambda_{i+1} + a\\y_{i+1} =& C \cdot x_{i+1} + D \cdot \lambda_{i+1} + F \cdot u(t_{i+1}) + e\end{aligned}\end{align} \]

It comes straightforwardly:

\[0 \leq y_{i+1} \, \perp \, \lambda_{i+1} \geq 0\]

These equations are summarized assuming that there is no source term and simplified relations as for the electrical oscillator with half-wave rectifier.

\[\begin{split}\begin{array}{ccc}
W & = & (I - h \theta A)^{-1} \\
x_{free} & = & W(I + h (1-\theta) A) \cdot x_{i} \\
x_{i+1} & = & x_{free} + h W r_{i+1} \\
r_{i+1} & = & B \cdot \lambda_{i+1} \\
y_{i+1} & = & C \cdot x_{i+1} + D \cdot \lambda_{i+1} \\
& 0 \leq y_{i+1} \, \perp \, \lambda_{i+1} \geq 0 &
\end{array}\end{split}\]

The integration algorithm with a fixed step is described here :

Algorthm 1: Integration of the electrical oscillator with half-wave rectifier through a fixed Moreau time stepping scheme

Require:\(R > 0 , L > 0 , C > 0\)

Require:Time parameters \(h,T,t_0\) and \(\theta\) for the integration

Require:Initial value of inductor voltage \(v_L = x_0(0)\)

Require:Optional, initial value of inductor current \(i_L = x_0(1)\) (default: 0)\(n_{step} = \frac{T - t_0}{h}\)

//

Dynamical system specification\(A = \left( \begin{array}{cc} 0 & \frac{-1}{C}\\\frac{1}{L} & 0 \end{array} \right)\)

//

Relation specification\(B = \left( \begin{array}{c}\frac{-1}{C}\\0\end{array} \right)\)

\(C = \left( \begin{array}{cc}-1 & 0\end{array} \right)\)

\(D = (R)\)

//

Construction of time independent operators

Require:\(I - h \theta A\) invertible \(W = (I - h \theta A)^{-1}\) \(M = D + h C W B\)//

Non-smooth dynamical system integration

for\(i=0\) to \(n_{step}-1\)do\[ \begin{align}\begin{aligned}x_{free} =& W (I + h (1 - \theta) A) x_i && \textrm{// Computation of $x_{free}$}\\q =& C \cdot x_{free} && \textrm{// Formalization of the one step LCP}\\(y_{i+1},\lambda_{i+1}) =& \textrm{solveLCP}(M,q) && \textrm{// One step LCP solving}\\x_{i+1} =& x_{free} + h W B \lambda_{i+1} && \textrm{// Computation of new state}\end{aligned}\end{align} \]

end for

We have used the SMASH simulator from Dolphin to perform a simulation of this circuit with a smooth model of the diode as given by Shockley’s law , with a classical one step solver (Newton-Raphson) and a choice between backward-Euler and trapezoidal integrators.

fig 3: Diodes characteristics from SPICE model depicts the static \(I(V)\) characteristic of two diodes with default SPICE parameters and two values for the emission coefficient \(N\): 1.0 (standard diode) and 0.25 (stiff diode).

The stiff diode is close to an ideal one with a threshold of 0.2 V.

fig 4: SMASH and SICONOS simulation results with backward Euler integration, 10 μs time step displays a comparison of the SMASH and SICONOS results with a backward Euler integration and a fixed time step of 10 μs. A stiff diode model was used in SMASH simulations. For fig 5: SMASH and SICONOS simulation results with trapezoidal integration, 10 μs time step, a trapezoidal integrator was used, yielding a better accuracy. One can notice that the results from both simulators are very close. The slight differences are due to the smooth model of the diode used by SMASH, and mainly to the threshold of around 0.2 V. Such a threshold yields small differences in the conduction state of the diode with respect to the ideal diode.