Going through each Step

Creating Transadmittance Matrix

It needs several steps to get the transadmittance matrices $[\tilde{Y}]$ and $[\hat{Y}]$ mentioned in equation (7.1). First the MNA matrix of the linear subcircuit (figure 7.1) is created (chapter 3.1) without the voltage sources $S_1$...$S_M$ and without the non-linear components. Note that all nodes must appear in the matrix, even those where only non-linear components are connected. Then the transimpedance matrix is derived by exciting one by one the port nodes of the MNA matrix with unity current. After that the transadmittance matrix is calculated by inverting the transimpedance matrix. Finally the matrices $[\tilde{Y}]$ and $[\hat{Y}]$ are filled with the corresponding elements of the overall transadmittance matrix.

Note: The MNA matrix of the linear subcircuit has $L$ nodes. Every node, that is connected to the non-linear subcircuit or/and is connected to a voltage source, is called "port" in the following text. So, there are $M+N$ ports. All these ports need to be differential ones, i.e. without ground reference. Otherwise problemes may occur due to singular matrices when calculating $[\tilde{Y}]$ or $[\hat{Y}]$.

Now this should be described in more detail: By use of the MNA matrix $[A]$, the $n$-th column of the transimpedance matrix $[Z]$ should be calculated. The voltage source at port $n$ is connected to node $i$ (positive terminal) and to node $j$ (negative terminal). This results in the following equation. (If port $n$ is referenced to ground, the -1 is simply omitted.)

$$[A]\cdot \begin{bmatrix}V_1\\ \vdots\\ V_L\\ \end{bmatrix} = \beg... ...\leftarrow i\text{-th row}\\ \\ \leftarrow j\text{-th row}\\ \\ \\ \end{matrix}$$ (7.3)

After having solved it, $Z_{1,n}$...$Z_{N+M,n}$ are obtained simply by subtraction of the node voltages:

$$Z_{m,n} = V_k - V_l$$ (7.4)

Here the voltage source at port $m$ is connected to node $k$ (positive terminal) and to node $l$ (negative terminal).

The next column of $[Z]$ is obtained by changing the right-hand side of equation (7.3) appropriately. As this has to be done $N+M$ times, it is strongly recommended to use LU decomposition.

As $[\tilde{Y}]$ is not square, problems encounter by trying to build its inverse matrix. Therefore, the following procedure is recommended:

• Create the transimpedance matrix for all ports (sources and interconnections).
• Compute the inverse matrix (transadmittance matrix).
• The upper left and upper right corner contains $[\tilde{Y}]$ and $[\hat{Y}]$.
• The lower left and lower right corner contains the transadmittance matrices to calculate the currents through the sources. They can be used to simplify the AC simulation at the very end.

One further thing must be mentioned: Because the non-linear components and the sources are missing in the linear MNA matrix, there are often components that are completely disconnected from the rest of the circuit. The resulting MNA matrix cannot be solved. To avoid this problem, shunt each port with a $100\Omega$ resistor, i.e. place a resistor in parallel to each non-linear component and to each source. The effect of these resistors can be easily removed by subtracting 10mS from the first main diagonal of the transadmittance matrix.

Starting Values

A difficult question is how to find appropriate start values for the harmonic balance simulation. It is recommended to first perform a DC analysis and start the algorithm with this result. In many situation (perhaps always) an even better starting point can be achieved by also using the result of a linear AC simulation. However with a large signal strength and strong non-linearities, convergence may still fail. Then, the following procedure might succeed: Perform HB simulation by applying half of the desired signal levels. If convergence is reached, the result can be used as start values for the simulation with the full signal levels. Otherwise the amplitude of the signals can be further decreased in order to repeat the above-mentioned procedure.

Solution algorithm

To perform a HB simulation, the multi-dimensional, non-linear function 7.2 has to be solved. One of the best possibilities to do so is the Newton method:

$$\textbf{V}_{n+1} = \textbf{V}_n - \textbf{J}_F (\textbf{V}_n)^{-1} \cdot \textbf{F} (\textbf{V}_n)$$ (7.5)

$$\Rightarrow \qquad \textbf{J}_F (\textbf{V}_n) \cdot \textbf{V}_{n+1} = \textbf{J}_F (\textbf{V}_n) \cdot \textbf{V}_n - \textbf{F} (\textbf{V}_n)$$ (7.6)

with $\textbf{J}_F$ being the Jacobian matrix. DC and transient simulation also use this technique, but here a problem appears: The derivatives of the component models are not given in frequency domain. Thus, the Jacobian must be calculated starting at the HB equation 7.2:

$$\boldsymbol{J}_F (\boldsymbol{V}_n) = \frac{d\boldsymbol{F} (\bol... ...N} + \boldsymbol{J}_{F,G} + j\cdot \boldsymbol{\Omega}\cdot\boldsymbol{J}_{F,Q}$$ (7.7)

So, two Jacobian matrices have to be built, one for the current $\boldsymbol{I}_G$ and one for the charge $\boldsymbol{Q}$. Both resulted from a Fourier Transformation. The two operations (Fourier Transformation and differentiation) are linear and thus, can be exchanged. Hence, the Jacobian matrices are built in time domain and transformed into frequency domain afterwards.

To obtain a practical algorithm of this procedure, the DFT is best written as matrix equation. By having a look at equation 15.180 and 15.181, it becomes clear how this works. The harmonic factors $\exp(j\omega_k t_n)$ build the matrix $\boldsymbol{\Gamma}$:

$$\qquad \boldsymbol{U}(j\omega) = \boldsymbol{\Gamma}\cdot \boldsymbol{u}(t)$$ (7.8)

$$\qquad \boldsymbol{u}(t) = \boldsymbol{\Gamma}^{-1}\cdot \boldsymbol{U}(j\omega)$$ (7.9)

with $\boldsymbol{u}$ and $\boldsymbol{U}$ being the vectors of the time and frequency values, respectively. Now, it is possible to transform the desired Jacobian matrix into frequency domain:

$$\boldsymbol{J}_{F,G} = \frac{\partial\boldsymbol{I}_G}{\partial\b... ...c{\partial\boldsymbol{i}}{\partial\boldsymbol{v}} \cdot\boldsymbol{\Gamma}^{-1}$$ (7.10)

Here $\boldsymbol{i}$ is a vector with length $K\cdot N$, i.e. first all time values of the first node are inserted, then all time values of the second node etc. The Jacobi matrix of $\boldsymbol{i}$ is defined as:

$$\boldsymbol{J}_{F,G}(\boldsymbol{u}) = \begin{bmatrix}\dfrac{\par... ...n}{\partial u_1} & \hdots & \dfrac{\partial i_n}{\partial u_n} \\ \end{bmatrix}$$ (7.11)

Hence this matrix consists of $K \times K$ blocks (one for each node) that are diagonal matrices with time values of the derivatives in it. (Components exists that create non-diagonal blocks, but these are so special ones that they do not appear in this document.)

The formula 7.10 seems to be quite clear, but it has to be pointed out how this works with FFT algorithm. With $\boldsymbol{\Gamma}^{-1} = (\boldsymbol{\Gamma}^{-1})^T$ (see equation 15.181) and $(\boldsymbol{A}\cdot\boldsymbol{B})^T = \boldsymbol{B}^T\cdot \boldsymbol{A}^T$, it follows:

$$\boldsymbol{J}_{F,G} = \boldsymbol{\Gamma}\cdot\frac{\partial\bol... ...\cdot \frac{\partial\boldsymbol{i}}{\partial\boldsymbol{v}} \right)^T \right)^T$$ (7.12)

So, there are two steps to perform an FFT-based transformation of the time domain Jacobian matrix into the frequency domain Jacobian:

1. Perform an FFT on every column of the Jacobian and build a new matrix $\boldsymbol{A}$ with this result, i.e. the first column of $\boldsymbol{A}$ is the FFTed first column of the Jacobian and so on.
2. Perform an IFFT on every row of the matrix $\boldsymbol{A}$ and build a new matrix $\boldsymbol{B}$ with this result, i.e. the first row of $\boldsymbol{B}$ is the IFFTed first row of $\boldsymbol{A}$ and so on.

As the Fourier transformation has to be applied to diagonal sub-matrices, the whole above-mentioned process can be performed by one single FFT. This is done by taking the $\partial\boldsymbol{i} / \partial\boldsymbol{v}$ values in a vector $\boldsymbol{J}_i$ and calculating:

$$\dfrac{1}{K}\cdot \left(\boldsymbol{J}_i\right)$$ (7.13)

The result is the first column of $\boldsymbol{J}_{F,G}$. The second column equals the first one rotated down by one element. The third column is the second one rotated down by one element etc. A matrix obeying this structure is called Toeplitz matrix.

So, finally the complete HB Newton iteration step can be written down. Putting 7.2 and 7.7 into 7.6 leads to

$$\textbf{J}_F (\textbf{V}_n) \cdot \textbf{V}_{n+1} = \textbf{J}_{... ...ol{\Omega}\cdot (\textbf{J}_{F,Q}\cdot\textbf{V}_n - \textbf{Q}) - \textbf{I}_S$$ (7.14)

This is important to notice, because many non-linear components cannot be processed at every bias point (see figure 3.5). These components create a new voltage estimate across their nodes, whereas the new estimated absolute voltages at their nodes are not known. Thus, the term $\textbf{J}_{F,G} \cdot \textbf{V}_n$ can only be created in one single step, leading to the vector $\textbf{I}_{G,J}$. Luckily, this procedure also saves computation time, as the matrix multiplication need not to be performed. The same is true for the term $\textbf{J}_{F,Q}\cdot\textbf{V}_n$, leading to the vector $\textbf{Q}_J$. So it is:

$$\textbf{J}_F \cdot \textbf{V}_{n+1} = \textbf{I}_{G,J} - \textbf{I}_G + j\cdot \boldsymbol{\Omega}\cdot (\textbf{Q}_J - \textbf{Q}) - \textbf{I}_S$$ (7.15)

Termination Criteria

Frequency components with very different magnitude appear in harmonic balance simulation. In order to detect when the solver has found an accurate solution, an absolute as well as relative criteria must be used on all nodes and at all frequencies. The analysis is regarded as finished if one of the criteria is satisfied.

The absolute and relative criteria write as follows:

$$\left\vert \tilde{I}_{n,k} + \hat{I}_{n,k} \right\vert < \varepsilon_{abs} \forall \quad n, k$$ (7.16)

$$2\cdot \left\vert \frac{\tilde{I}_{n,k} + \hat{I}_{n,k}} {\tilde{I}_{n,k} - \hat{I}_{n,k}} \right\vert < \varepsilon_{rel} \forall \quad n, k$$ (7.17)

where $\tilde{I}_{n,k}$ is the current of the linear circuit partition for node $n$ and frequency $k$ and $\hat{I}_{n,k}$ is the current of the non-linear circuit partition.