# Model-Free Closed-Loop Stability Analysis: A Linear Functional Approach

Adam Cooman, Fabien Seyfert, Martine Olivi, Sylvain Chevillard, Laurent Baratchart

► To cite this version:

Adam Cooman, Fabien Seyfert, Martine Olivi, Sylvain Chevillard, Laurent Baratchart. Model-Free Closed-Loop Stability Analysis: A Linear Functional Approach. IEEE Transactions on Microwave Theory and Techniques, 2018, 66 (1), pp.73-80. 10.1109/TMTT.2017.2749222. hal-01381731

HAL Id: hal-01381731

https://inria.hal.science/hal-01381731v1

Submitted on 11 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Model-Free Closed-Loop Stability Analysis:
A Linear Functional Approach

Adam Cooman, Fabien Seyfert, Martine Olivi, Sylvain Chevillard and Laurent Baratchart*

Abstract

Performing a stability analysis during the design of any electronic circuit is critical to guarantee
its correct operation. A closed-loop stability analysis can be performed by analysing the impedance
presented by the circuit at a well-chosen node without internal access to the simulator. If any of the
poles of this impedance lie in the complex right half-plane, the circuit is unstable. The classic way
to detect unstable poles is to fit a rational model on the impedance.

In this paper, a projection-based method is proposed which splits the impedance into a stable and
an unstable part by projecting on an orthogonal basis of stable and unstable functions. When the
unstable part lies significantly above the interpolation error of the method, the circuit is considered
unstable. Working with a projection provides one, at small cost, with a first appraisal of the unstable
part of the system.

Both small-signal and large-signal stability analysis can be performed with this projection-based
method. In the small-signal case, a low-order rational approximation can be fitted on the unstable
part to find the location of the unstable poles.

Frequency domain simulation methods, like Harmonic Balance (HB) or a DC analysis, impose a
structure on the obtained solution of the circuit [21]: The DC analysis only allows for a fixed solution,
while HB imposes a frequency grid. Any circuit solution that requires more than the imposed frequencies,
e.g. an extra oscillation not on the imposed grid, cannot be represented in the constrained frequency grids
of DC and HB. The simulator will still find a valid solution, but the obtained orbit will be locally unstable:
it cannot recover from small perturbations and will be physically unobservable in the circuit [22]. It is
therefore necessary to perform a local stability analysis on each of the circuit solutions obtained with a
DC and HB analysis [21].

Over the years, several methods have been developed to determine the local stability of a circuit
solution. Some techniques, like the analysis of the characteristic system [22], require access to the
simulator. Open-loop techniques, like the analysis of the normalised determinant [21], require access
to the intrinsic device models. These classic techniques are therefore hard to implement in commercial
simulators.

Closed-loop stability analysis methods can easily be applied as a post-processing step without any
internal knowledge of the circuit and can be used in commercial simulators. This is the reason why
they have attracted a large interest lately [5, 10, 21]. A closed-loop local stability analysis analyses the
linearisation of the circuit around the orbit to determine the stability of that orbit. If the linearised
circuit has poles in the complex right half-plane, the orbit is unstable.

The poles of the linearisation around the circuit orbit cannot be obtained directly. Instead a Frequency
Response Function (FRF) of the linearised circuit is obtained with small-signal simulations on a discrete
set of frequencies. The closed-loop stability analysis then verifies whether the underlying FRF has any
poles in the complex right half-plane.

In a pole-zero stability analysis, a rational approximation is fitted on the FRFs. If the rational approximation contains poles in the complex right half-plane, the solution is unstable.

The FRF of circuits with distributed elements, like transmission lines, cannot be fitted exactly by
a rational function over a wide frequency range. A good approximation to the FRF can be obtained

<footnote>*Adam Cooman, Fabien Seyfert, Martine Olivi, Sylvain Chevillard and Laurent Baratchart are with APICS at INRIA, Sophia Antipolis. e-mail: adam.cooman@inria.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: A Small-signal current source is connected to a well-chosen node in the circuit under test to perform the local stability analysis.</figcaption></figure>

using modern rational approximation techniques and rational functions with a very high order, but the obtained rational function will contain spurious poles [20] which could lead to false positives in a stability analysis [1].

A heuristic approach to overcome the issues introduced by distributed elements, proposed in [1], is to cut the frequency band into smaller pieces and to use low-order local rational approximations to determine the stability of each part of the total FRF. This algorithm is commercially available in the STAN tool [7, 13] and it was successfully applied to several examples [2, 6, 17, 23].

In this paper, we propose an alternative closed-loop stability analysis method. In the proposed stability analysis, the FRF is split into a stable and unstable part by projection onto the orthogonal basis of stable and unstable functions. If a significant part of the FRF is projected onto the unstable basis functions, the circuit solution is unstable. Calculating the projection boils down to calculating a Fourier transform when the FRF is transformed from the complex plane to the unit disc. Using the Fast Fourier Transform (FFT), this can be done fast and in a numerically robust way.

Functional projection onto a stable and unstable basis is a linear operation, simple to implement, and no optimisation step is required. No model-order or maximum approximation error needs to be specified. The parameters in the projection method are the frequency range on which the FRF is determined and the amount of simulation points. When the amount of simulated points is too low, an interpolation error is present in the result of the method. It is shown that the level of this error can easily be estimated and used to correctly choose the amount of needed simulation points.

Once the unstable part of the FRF has been obtained, it is compared to the level of the interpolation error to determine whether the unstable part is significant or not. This final step can be done visually, or a significance threshold can be chosen by the user, both will require some experience with the method.

A final benefit of the projection-based approach is that it may help exploiting the fact that the unstable part is rational in a small-signal stability analysis [4]. The unstable part can therefore be approximated by a rational function without influence of the distributed elements, which are projected onto the stable part of the FRF.

The following of the paper is structured as follows: First, the simulation set-up used to determine the FRF of the linearised circuit is discussed (Section 1). Then, the details of the functional projection are provided (Section 2). In Section 3, the method is applied to four examples: First, an artificial example is considered. Then, the small-signal stability of two amplifiers is investigated and finally, the method is applied to investigate the large-signal stability of a circuit.

# 1 Determining the frequency responses

In this paper, the (trans)impedance presented by the circuit to a small-signal current source will be used as FRF (Fig. 1). In the remainder of this paper, it will be assumed that the unstable poles are observable in the FRF. To reduce the chance of missing an instability in the circuit due to a pole-zero cancellation, many different FRFs can be analysed one-by-one. Having a fast method to determine stability of a single FRF is therefore critical to a robust stability analysis.

The FRF of the linearised circuit is obtained by first placing the circuit in the required orbit, using either a DC or HB analysis and running a small-signal simulation around this orbit.

In a small-signal stability analysis, the stability of the DC solution of the circuit is investigated, so the FRF of the linearised circuit is obtained with an AC simulation. The impedance of the circuit is
then obtained as:

$$Z_{mn}(j\omega) = \frac{V_m(j\omega)}{I_n(j\omega)} \qquad (1)$$

where $I_n(j\omega)$ is the small-signal current injected into the selected node n and $V_m(j\omega)$ is the voltage response of the circuit measured at node m in the circuit.

In a large-signal stability analysis, the stability of a large-signal solution of the circuit is investigated. The circuit is driven by a periodic continuous-wave excitation at a pulsation $\omega_{ex}$ and the circuit solution is obtained with a HB simulation.

The FRF of the linearised system around the HB orbit is obtained with a mixer-like simulation<fnref n="1" />. As the small-signal will mix with the large signal, several transfer impedances with a different frequency translation are obtained:

$$Z_{mn}^{[b]}(j\omega) = \frac{V_m(j\omega + b j\omega_{ex})}{I_n(j\omega)} \quad b \in \mathbb{Z}$$

The stability analysis now needs to determine whether the obtained impedances have poles in the right half-plane. The stability analysis of a large-signal orbit doesn't differ much from the analysis of a DC solution [5]. The small-signal stability analysis can be considered a special case where only $Z_{mn}^{[0]}(j\omega)$ is analysed.

## 2 Stable/unstable projection

The projection described here has been used before to perform stable interpolation and extrapolation of FRF data [16]. With slight modifications, it can be turned into a full-blown stability analysis. We first start with a brief introduction to the notion of Hardy spaces.

The Hardy space $H^2(\mathbb{C}^+)$, is defined as the set of all functions $g$ defined on $\mathbb{C}^+$ such that:

* $\forall z \in \mathbb{C}^+$, $g$ is holomorphic at $z$

* $\sup_{x>0} \int_{-\infty}^{+\infty} |g(x + j\omega)|^2 d\omega < \infty$

A classical result [9, 18] states that every function $g \in H^2$ admits a limiting function $G(j\omega)$ defined on the imaginary axis $j\mathbb{R}$. The latter is obtained by taking the limit of $g(z)$ when $z$ tends non tangentially toward $j\omega$. Moreover $\forall z \in \mathbb{C}^+$, $g(z)$ is equal to the poisson integral of $G$, that is:

$$g(z = x + jy) = \int_{-\infty}^{+\infty} G(j\omega) \frac{x}{x^2 + (y - \omega)^2} d\omega. \qquad (2)$$

The holomorphic nature of $g$ ensures that it is also the Cauchy integral of its boundary value $G$:

$$g(z = x + jy) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} G(j\omega) \frac{1}{j\omega - z} d\omega. \qquad (3)$$

There is a one to one linear correspondence between the Hardy functions $g$ and its boundary value function $G$. Using this identification, $H^2(\mathbb{C}^+)$ becomes a subspace of the Hilbert space $L^2(j\mathbb{R})$ of square integrable functions on $j\mathbb{R}$. A direct consequence of (2) is that $H^2(\mathbb{C}^+)$ is closed in $L^2(j\mathbb{R})$ and therefore admits an orthogonal complement. An important result [9] asserts that, $(H^2(\mathbb{C}^+))^\perp = H^2(\mathbb{C}^-)$ where $H^2(\mathbb{C}^-)$ is defined exactly as $H^2(\mathbb{C}^+)$ above by replacing $\mathbb{C}^+$ by $\mathbb{C}^-$ and taking the supremum over $x < 0$. We therefore have that

$$L^2(j\mathbb{R}) = H^2(\mathbb{C}^+) \oplus H^2(\mathbb{C}^-)$$

This decomposition asserts that any square integrable function on $j\mathbb{R}$ decomposes uniquely as the sum of the traces on the imaginary axis, of an analytic function in the right half-plane and a function analytic in the left half-plane. The projection on $H^2(\mathbb{C}^+)$ defines the stable part of the function. The projection

<footnote n="1">¹In Keysight's Advanced Design System (ADS), this mixer-like simulation is called a Large-Signal Small-Signal (LSSS) analysis.</footnote>
onto $H^2(\mathbb{C}^-)$ is the unstable part. As an example, consider $P/Q$ a strictly proper ($\deg(P) < \deg(Q)$) rational function devoid of poles on the imaginary axis. We write its partial fraction expansion as,

$$ \frac{P(s)}{Q(s)} = \sum_{i \in I^+} \sum_{k=1}^{k_i} \frac{a_{i,k}}{(s - \lambda_i)^k} + \sum_{i \in I^-} \sum_{k=1}^{k_i} \frac{a_{i,k}}{(s - \lambda_i)^k} $$

where the $\lambda'_i$s with $i \in I^-$ are poles belonging to $\mathbb{C}^-$, and the ones with $i \in I^+$ belong to $\mathbb{C}^+$. The strict properness of $P/Q$ ensures its square integrability on $j\mathbb{R}$. By unicity its stable part obtained after projection on $H^2(\mathbb{C}^+)$ is found to be,

$$ \sum_{i \in I^-} \sum_{k=1}^{k_i} \frac{a_{i,k}}{(s - \lambda_i)^k}, $$

while its unstable part is

$$ \sum_{i \in I^+} \sum_{k=1}^{k_i} \frac{a_{i,k}}{(s - \lambda_i)^k}. $$

In the general case the projection boils down to calculating certain inner products of $Z_{mn}^{[b]}(j\omega)$ with the basis functions $B_k$, which form an orthogonal basis of $L^2(j\mathbb{R})$

$$ c_k = \langle Z_{mn}^{[b]}(j\omega), B_k \rangle = \int_{-\infty}^{\infty} Z_{mn}^{[b]}(j\omega) \overline{B_k(j\omega)} d\omega \quad (4) $$

$$ B_k(s) = - \sqrt{\frac{\alpha}{\pi}} \frac{(s - \alpha)^k}{(s + \alpha)^{k+1}} \quad k \in \mathbb{Z} \quad (5) $$

The overbar $\bar{\cdot}$ indicates the complex conjugate. $\alpha$ is a positive constant used for scaling. All $B_k$ with $k \ge 0$ create a basis for the stable part, while the $B_k$ with negative $k$ form a basis for $H^2(\mathbb{C}^-)$. Once the $c_k$ coefficients are calculated, the stable and unstable parts are easily recovered by calculating

$$ Z_{\text{stable}}(j\omega) = \sum_{k=0}^{\infty} c_k B_k(j\omega) \quad (6) $$

$$ Z_{\text{unstable}}(j\omega) = \sum_{k=1}^{\infty} c_{-k} B_{-k}(j\omega) \quad (7) $$

The inner product in (4) runs over all frequencies while the impedance function $Z_{mn}^{[b]}(j\omega)$ is only known over a frequency range $f = [f_{\min}, f_{\max}]$. To impose the finite frequency band on the data, the impedance is filtered before the analysis

$$ Z_f(j\omega) = Z_{mn}^{[b]}(j\omega) H(j\omega) \quad (8) $$

The filter $H(j\omega)$ is a high-order elliptic lowpass filter with its first transmission zero placed at $f_{\max}$ that imposes band limitation. This filter will stabilise poles close to its cutoff frequency, so $f_{\max}$ should be chosen well beyond the maximum frequency at which the circuit can become unstable. $f_{\min}$ should be placed very close to DC. If $f_{\min}$ can't be close to DC, a bandpass filter should be used for $H(j\omega)$. The filtering ensures that $Z_f(j\omega) \in L^2(j\mathbb{R})$ by suppressing anything outside of the frequency band of interest. The smooth decay to zero of $Z_f(j\omega)$ at the edges of the frequency interval will avoid instabilities to pop up due to the discontinuity of $Z_{mn}^{[b]}(j\omega)$.

In this paper, we use an elliptic filter of order 10 to filter the data. The filter has one transmission zero which is placed exactly at $f_{\max}$.

## 2.1 Transforming to the unit circle

Working with the basis functions defined in (5) is troublesome from a numeric point of view. Performing the projection when working on unit disc yields better results [16]. The mapping from the complex plane to the unit circle is performed with the Möbius mapping visualised in Fig. 2:

$$ Z_f(s) \underset{\text{Möb}}{\mapsto} Z_f^{\text{disc}} = \sqrt{\pi\alpha} \frac{2}{z-1} Z_f \left( \alpha \frac{1+z}{1-z} \right) \quad (9) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2: The Möbius transform used in the analysis maps the complex right-half plane to the inside of the unit disc. DC is mapped to -1. When $\alpha = 2\pi f_{max}(\sqrt{2}-1)$, $f_{max}$ is mapped to $e^{j\frac{\pi}{4}}$.</figcaption></figure>

Our mapping of choice converts square integrable functions on the frequency axis into square integrable functions of the same norm on the unit circle. Square integrable functions which are analytic on the right half-plane are mapped onto square integrable functions which are analytic inside the unit disc. Appendix 4 shows that the basis functions $B_k$ map onto powers of z

$$B_k(s) \stackrel{\text{Möb}}{\mapsto} B_k^{\text{disc}}(z) = z^k \qquad (10)$$

Projecting on this basis boils down to calculating the Fourier series of $Z_f^{\text{disc}}(z)$, with coefficients given by

$$c_k = \frac{1}{2\pi} \int_{0}^{2\pi} Z_f^{\text{disc}}(e^{j\theta}) e^{-jk\theta} d\theta \qquad (11)$$

This Fourier series can be calculated in a numerically efficient way using the Fast Fourier Transform (FFT). The FFT requires that the $\theta$-values are linearly spaced between 0 and $2\pi$. Due to the mapping from the complex plane to the unit disc, the samples will not satisfy this constraint. A simple interpolation, can be used to obtain $Z_f^{\text{disc}}$ on the $\theta$-values required to perform the FFT.

The interpolation can introduce artefacts in the unstable part if the FRF is not sampled on a sufficiently dense frequency grid, which will be shown on an example later. The level of this interpolation error can be estimated by interpolating the FRF using only the data at the even data points. The difference between the original and the interpolated data for the odd data points will give an indication of the interpolation error encountered in the stability analysis. When the unstable part lies significantly above the interpolation error, we can conclude that the original impedance is unstable.

The threshold to determine when the unstable part lies significantly above the level of the interpolation error will depend on the simulation set-up for the circuit. In the examples that follow, a level of 20 dB has been used as a threshold. A more strict threshold could be used at the cost of requiring a more dense frequency grid and an increased simulation time. When measured components are used in the simulation set-up, the noise level in those measurements should be taken into account to choose the correct threshold.

## 2.2 Summary of the projection method

The stable and unstable parts of a FRF are determined using the following steps:

1. Multiply the FRF with a high-order filter as in (8)

2. Transform the filtered FRF to the unit disc using (9)

3. Interpolate the transformed FRF to a linear grid and use the FFT to calculate the $c_k$ coefficients

4. Reconstruct the stable and unstable part using (6-7)

## 2.3 Obtaining the unstable poles

A function is meromorphic on $\mathbb{C}$ if it is holomorphic on $\mathbb{C}$ but on a countable number of isolated poles. The function $\tanh(\omega)$ is for example meromorphic, having infinitely many isolated poles on the imaginary axis placed at $j\pi/2 + jk\pi$ ($\forall k \in \mathbb{Z}$) and being analytic elsewhere. In a small-signal stability analysis of a
<figure><img src="image_4.png" /><figcaption>Fig. 3: Pole-zero map of the test system. There are 202 poles (+) and 200 zeroes (o). The two poles placed in the right half-plane are easily recovered after the projection by fitting a low-order model on the unstable part (□).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4: The stable/unstable projection splits the original frequency response (-) in a stable part indicated with (-) and an unstable part indicated with (–)</figcaption></figure>
circuit composed of lumped elements, transmission lines and active devices modelled by negative resistors, the impedances can be shown to be meromorphic functions of the frequency. Under the additional realistic assumption that active elements can only deliver power over a finite bandwidth, the impedances are proven to possess only finitely many unstable poles in $C^+$ [4]. Under the generic condition that $Z_{mn}^{[b]}(j\omega)$ is devoid of poles on the imaginary axis, and that the filtering function $H(j\omega)$ decays strongly enough in order to render $Z_{mn}^{[b]}(j\omega)H(j\omega)$ square integrable, we conclude that the unstable part of $Z_{mn}^{[b]}(j\omega)H(j\omega)$ is a rational function. Its poles coincide with the unstable poles of $Z_{mn}^{[b]}(j\omega)$: note here that the multiplication by $H(j\omega)$ does not add any unstable pole as $H(j\omega)$ is stable. This means that most of the complexity of the frequency response, like the delay, will be projected onto the stable part, while the unstable part can easily be approximated by a low-order rational model to recover the unstable poles<fnref n="2" />.

Classic rational approximation tools can be used to approximate the unstable part and determine the unstable poles. When multiple frequency responses are analysed simultaneously, an approximation method suited for approximation of rational matrices is preferred [15]. In our current implementation, Kung's method [11, 14] is used to estimate the poles of the unstable part of a single FRF at a time. Alternatively, more sophisticated rational approximation engines like RARL2 [15] can be used to recover and track unstable poles.

Compared to working with a high-order rational approximation of the total impedance of the circuit, the 'split-first, approximate later' approach proposed here could be a faster and easier method to recover the unstable poles. The post-processing will be faster, but the amount of points required to obtain a sufficiently low interpolation error might be higher than the amount of points required for a tool based on rational approximation.

When the circuit is stable, designers often require information about critical stable poles, to determine how far the circuit is from instability and to track the location of the poles as the circuit varies. In its current form, the projection-based analysis does not simplify finding the location of the stable poles. To perform such an analysis, methods based on local modelling may still be required.

## 3 Examples

The stability analysis will now be applied to four different examples. The first is an artificial example generated in Matlab on which we can demonstrate that the unstable poles in the circuit are recovered perfectly. In the second example, an unstable balanced amplifier is analysed to show that the method works for RF circuits. The third example is a two-stage GaN Power Amplifier (PA). In the final example, a large-signal stability analysis is performed to verify the stability of a circuit orbit obtained in a HB simulation.

All simulations were performed in Keysight's Advanced Design System (ADS) and the post-processing was performed in Matlab.

### 3.1 Example 1: Random state space system

As a first example, the stability analysis is applied to a random system of order 202 generated with the `rss` function from Matlab<fnref n="3" />. The test system has an unstable pole pair at 1 GHz, as can be seen on its pole-zero map (Fig. 3). A zero is placed close to the unstable poles. This makes that the unstable poles are difficult to observe in the FRF. To introduce delay in the test system, a time delay of 2 ns is added to the system. The frequency response of the system is calculated on 5000 linearly spaced frequency points between 0 Hz and 5 GHz and is shown in green in Fig. 4. The obtained stable and unstable parts after projection are shown in red and blue on the same figure. The maximum interpolation error is very low in this example (−120 dB). We will focus on the effect of the interpolation error in more detail in example 2.

The obtained unstable part peaks at 1 GHz, which matches the location of the unstable pole pair of the system. Note also that the obtained $Z_{\text{unstable}}$ is very simple: it is clearly a second-order system.

<footnote n="2">²Interpolation error will be present in the obtained unstable part, but its influence can be minimised by weighting the rational estimator with the obtained interpolation error level.</footnote>

<footnote n="3">³The `rss` function in Matlab returns models with poles and zeroes around 1Hz. The example here was scaled up in frequency to represent an RF circuit.</footnote>
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5: Simulation set-up and photograph of the balanced amplifier. In the simulations for the stability analysis, the amplifier is excited at the collector of one of its transistors. All transmission lines in the circuit are 50Ω lines unless stated otherwise. The length of the transmission lines is given in millimetres. The TLINP model was used for the transmission lines with $\epsilon_r = 6.15$, $\tan(\delta) = 0.003$ and conductor losses $A = 2.5$ dB/m.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6: The separation of the impedance (–) into the stable part (–) and unstable part (–) reveals the instability around 1.46 GHz. The interpolation error is shown with (–). In the plot on the left, there are artefacts present in the unstable part due to interpolation of the coarsely obtained impedance data. When the impedance of the balanced amplifier is simulated on a finer frequency grid, the artefacts disappear (right).</figcaption></figure>

Most of the complexity of the frequency response, including the delay, is projected onto the stable part. This observation supports the proposed approach of estimating a rational model only after projection. A good fit was obtained with a rational model that consisted of two unstable poles and a single zero. The two poles obtained with a rational approximation of $Z_{\text{unstable}}$ coincide exactly with the unstable poles in the circuit as is shown in Fig. 3.

## 3.2 Example 2: Balanced amplifier

The second example is a balanced amplifier built as a student project for operation around 3.4 GHz (Fig. 5). Two BFP520 transistors were used to construct the amplifier. During measurement, the design oscillated around 1.43 GHz when terminated with 50 Ω, so the circuit is a good candidate to verify the proposed method to find the instability in simulations.

The small-signal current source was connected to the collector of the top transistor. The BFP520 has a $f_T$ of 45 GHz, so the maximum frequency for the simulation was set to 50 GHz. The impedance of the circuit was determined starting from DC in 10 MHz steps. The obtained impedance is shown in green in Fig. 6, the obtained stable and unstable parts are also shown in the same figure. The instability around 1.46 GHz is detected, but also some artefacts can be observed in the obtained unstable part at higher frequencies. The high level of the interpolation error at the frequency of these artefacts indicates that they are due to the interpolation in step 3 of the stable/unstable projection. At the frequency of the detected instability, the unstable part lies about 30 dB above the error level, which indicates that the instability is not an interpolation artefact.
To confirm that the artefacts are caused by the interpolation, a second simulation was run, but now 1 MHz steps were used instead of 10 MHz. The stable/unstable projection of the denser frequency response data (Fig. 6) still predicts the instability around 1.46 GHz, but the artefacts in the unstable part at higher frequencies are gone. The maximum of the interpolation error went down to -80 dB(Ω).

### 3.3 Example 3: Two-stage Power Amplifier

<figure><img src="image_9.png" /><figcaption>Fig. 7: Microphotograph of the MMIC. The stabilisation resistor $R_{stab}$ is indicated in red.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8: Impedance seen at the gate of the first stage of the PA for $R_{stab} = 500Ω$. Its obtained stable and unstable parts clearly indicate that the DC solution of the amplifier is unstable. The interpolation error (-) is quite high due to the low amount of simulation points.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9: Impedance presented by the PA at the gate of the transistor in the first stage for $R_{stab} = 25Ω$. The obtained stable and unstable parts indicate that the DC solution is now stable. The interpolation error is shown with (-).</figcaption></figure>

As a third example, we consider the small-signal stability analysis of an X-band PA designed in the 0.25µm GaN HEMT technology GH25-10 of UMS [8]. The circuit and its design are described in great detail in [23]. The resulting MMIC is shown in Fig. 7.

The PA is a two-stage design where the second stage consists of two branches with each two transistors in parallel. In simulation, the second stage of the PA demonstrated an odd-mode instability [23], so a stabilisation resistor was added between the drains of the top and bottom halves of the second stage of the PA (as indicated in the Fig.).
The simulation of the complete PA was performed in ADS. The passive structures in the circuit were simulated with EM simulations in Momentum and combined with the non-linear transistor models afterwards. To verify the stability of the amplifier, the circuit impedance was determined at the gate of the top most transistor of the second stage.

The obtained impedance for $R_{\text{stab}} = 500\Omega$ is shown in Fig. 8. The impedance is simulated on 945 logarithmically spaced points between 1MHz and 50GHz. Because 1MHz is not sufficiently close to DC, a bandpass filter was used in the stability analysis. Due to the low amount of data points in the resonances of the circuit, a Padé interpolation was used in the stability analysis. The resulting stable and unstable parts are shown in blue and red on the same figure. It is clear that the circuit is unstable for $R_{\text{stab}} = 500\Omega$. The unstable part peaks around 9.5GHz and lies about 40 dB above the interpolation error level.

The odd-mode instability can be resolved by decreasing the resistance of $R_{\text{stab}}$ [23]. In a second stability analysis, we determined the stability of the PA for $R_{\text{stab}} = 25\Omega$. The results of this second analysis are shown in Fig. 9. The obtained unstable part coincides with the level of the interpolation error, which indicates that the circuit is now stable.

## 3.4 Example 4: R-L-diode circuit

<figure><img src="image_12.png" /><figcaption>Fig. 10: The circuit R-L-diode circuit is excited with a small-signal current source at the diode.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11: The bifurcation diagram of the R-L-diode circuit shows that a period-doubling occurs for input amplitudes higher than 0.8 V.</figcaption></figure>

The final example in this paper shows that the stability analysis can also be used to determine the stability of HB simulations of the R-L-diode circuit shown in Fig. 10. The circuit is based on [3], but a realistic diode model was used to represent the diode in the circuit instead of the three equations provided in the original paper.

The circuit is excited by a single-tone voltage source with an amplitude $V_{\text{in}}$ and a frequency of 100 kHz. Because the diode has a transit-time of 4 µs, the circuit generates period-doubling solutions starting from sufficiently high amplitudes $V_{\text{in}}$. For even higher $V_{\text{in}}$, the circuit will create chaotic solutions.

To visualise this behaviour, a bifurcation diagram is constructed using time-domain simulations in the same way as is described in [3]: For every value of $V_{\text{in}}$, 1030 periods of 100 kHz are simulated and the final 30 periods are sampled every $1/100$ kHz. If the circuit solution is periodic with the same period as the input source, all 30 sampled points will fall on top of each-other. If a period-doubling occurs in the circuit, two different values will be obtained.
<figure><img src="image_14.png" /><figcaption>Fig. 12: The results of the stability analysis of $Z_{mn}^{[-1]}(j\omega)$, $Z_{mn}^{[0]}(j\omega)$ and $Z_{mn}^{[+1]}(j\omega)$ in the two HB simulations show that, for $V_{in} = 0.5$ V, the orbit is stable, but for $V_{in} = 1.5$ V, the orbit is unstable. The interpolation error is shown with (-).</figcaption></figure>

The obtained bifurcation diagram for our R-L-diode example is shown in Fig. 11. It is clear that a period-doubling occurs for $V_{in}$ higher than 0.8 V. Starting from 1.8 V the period quadruples. For the highest input amplitudes, a chaotic solution is obtained.

If this R-L-diode circuit is simulated with HB, the circuit solution is constrained to harmonics of 100 kHz. For input amplitudes higher than 0.8 V, where the circuit wants to go to a period-doubling solution, the constrained HB solution will be locally unstable.

We run two HB simulations on this circuit. Both HB simulations have a base frequency of 100 kHz and an order of 10. In the first simulation, $V_{in}$ is set to 0.5 V, which will result in a stable orbit. The second simulation has a $V_{in}$ of 1.5 V, which will cause the orbit to be unstable.

The frequency response of the circuit around the HB solution is obtained with a mixer-like simulation, as explained in the introduction of this paper. The small-signal excitation was swept in both cases on a linear frequency grid starting from (1 kHz + 1 Hz) up to (2 MHz + 1 Hz) in 1 kHz steps. The 1 Hz was added to the start and stop values of the sweep to avoid overlap with the tones of the HB simulation.

The mixer-like simulation in ADS uses Single-Sideband (SSB) current excitations $i(t) = e^{j\omega t}$, which causes the obtained frequency responses $Z_{mn}^{[b]}(j\omega)$ with $b \neq 0$ to be non-Hermitian:

$$Z_{mn}^{[b]}(j\omega) \neq \overline{Z_{mn}^{[b]}(-j\omega)}$$

An alternative representation can make $Z_{mn}^{[b]}(j\omega)$ Hermitian by transferring to a sine and cosine basis from the exponential basis [12, 19]

$$Z_{mn}^{[b]}(j\omega) = \frac{1}{2} [Z_{mn}^{[b]}(j\omega) + \overline{Z_{mn}^{[-b]}(j\omega)}] \\ Z_{mn}^{[-b]}(j\omega) = \frac{j}{2} [Z_{mn}^{[b]}(j\omega) - \overline{Z_{mn}^{[-b]}(j\omega)}]$$

$Z_{mn}^{[-1]}(j\omega)$, $Z_{mn}^{[0]}(j\omega)$ and $Z_{mn}^{[+1]}(j\omega)$ are then analysed with the stable/unstable projection method. The results are shown in Fig. 12. The HB solution obtained for $V_{in} = 0.5$ V is clearly stable: its unstable part is more than 70 dB smaller than its stable part.

In the case for $V_{in} = 1.5$ V, the solution is clearly unstable as the unstable part lies far above the stable part of the frequency response. Note that the lowest-frequency peak in the unstable part is located around 50 kHz and that copies of the resonance are found at 150 kHz, 250 kHz,... This behaviour is to be expected and indicates that the circuit wants to go to a period-doubling solution.

During the stability analysis of a periodic orbit, the unstable part will contain both the unstable base pole and all its higher-order copies. The unstable part will be simple, just like in the small-signal case, and it will be possible to approximate it by a finite set of base poles. Due to the infinite amount of higher-order copies however, it will not be possible to approximate it by a low-order rational approximation as is the case in the stability analysis of a DC solution.
# 4 Conclusion

This paper introduces a closed-loop local stability analysis without using a rational approximation. Instead, the impedance functions are split into a stable and unstable part by projecting onto an orthogonal basis. Transforming the problem to the unit disc allows to calculate this projection with the FFT which makes the projection-based stability analysis very fast. In a small-signal stability analysis, once the unstable part is obtained, a low-order rational model can be used to find the unstable poles in the circuit.

Due to the model-free nature of the proposed method, it is a very simple method to use: no choice of model order or approximation error needs to be made. The only requirements of the projection-based stability analysis are that the frequency responses are sampled on a sufficiently dense frequency grid and that the maximum frequency of the simulations is large enough. When the circuit impedance is simulated on a too coarse frequency grid, a large interpolation error is introduced in the results. The level of this interpolation error can easily be determined and used to improve the accuracy of the method.

Once the stable and unstable parts of the impedance are obtained with a sufficiently low interpolation error, the obtained unstable part can be compared to the interpolation error to determine whether it is significant or not. From experience, we found that, when the unstable part lies more than 20 dB above the interpolation error level, the circuit can be considered unstable. Further work is to be done towards automated decision making regarding stability.

The stable/unstable projection has been successfully applied to both the stability analysis of DC and large-signal solutions of RF circuits.

# Appendix: mapping of the basis functions onto the unit disc

Applying transform (9) to the basis functions of the complex plane (5) yields the following:

$$B_k^{\text{disc}}(z) = \sqrt{\pi\alpha} \frac{2}{z-1} B_k\left(\alpha \frac{1+z}{1-z}\right) = -\sqrt{\pi\alpha} \frac{2}{z-1} \sqrt{\frac{\alpha}{\pi} \left(\alpha \frac{1+z}{1-z} - \alpha\right)^k} = z^k$$

# Acknowledgement

This research was partly supported by the French space agency CNES and partly by the Flemish Agency for Innovation by Science and Technology (IWT-Vlaanderen). We are also thankful to Juan-Marie Collantes (UPV) for fruitful discussions on the topic of closed loop stability analysis. We would like to thank Kurt Homan, Johan Nguyen and Dries Peumans for the design and measurement of the balanced amplifier. Finally, we would like to thank Marc van Heijningen for providing the data of the MMIC PA.

# References

[1] A. Anakabe, N. Aylloon, J.M. Collantes, A. Mallet, G. Soubercaze-Pun, and K. Narendra. Automatic pole-zero identification for multivariable large-signal stability analysis of rf and microwave circuits. In *Microwave Conference (EuMC), 2010 European*, pages 477–480, Sept 2010.

[2] N. Ayllon, J. M. Collantes, A. Anakabe, I. Lizarraga, G. Soubercaze-Pun, and S. Forestier. Systematic approach to the stabilization of multitransistor circuits. *IEEE Tran. on Microwave Theory and Techniques*, 59(8):2073–2082, 2011.

[3] A. Azzouz, R. Duhr, and M. Hasler. Transition to chaos in a simple nonlinear circuit driven by a sinusoidal voltage source. *IEEE Transactions on Circuits and Systems*, 30(12):913–914, Dec 1983.

[4] L. Baratchart, S. Chevillard, and F. Seyfert. On transfer functions realizable with active electronic components. hal-01098616, 2014. [Research Report] RR-8659.

[5] J.M. Collantes, I. Lizarraga, A. Anakabe, and J. Jugo. Stability verification of microwave circuits through floquet multiplier analysis. In *Circuits and Systems, 2004. Proceedings. The 2004 IEEE Asia-Pacific Conference on*, volume 2, pages 997–1000 vol.2, Dec 2004.
[6] J.M. Collantes, N. Otegi, A. Anakabe, N. Ayllon, A. Mallet, and G. Soubercaze-Pun. Monte-carlo stability analysis of microwave amplifiers. In *Wireless and Microwave Technology Conference (WAMICON), 2011 IEEE 12th Annual*, pages 1–6, April 2011.

[7] AMCAD engineering. *STAN tool - Stability Analysis*.

[8] D. Floriot, H. Blanck, D. Bouw, F. Bourgeois, M. Camiade, L. Favede, M. Hosch, H. Jung, B. Lambert, A. Nguyen, K. Riepe, J. Splettstosse, H. Stieglauer, J. Thorpe, and U. Meiners. New qualified industrial algan/gan hemt process: Power performances amp; reliability figures of merit. In *2012 7th European Microwave Integrated Circuit Conference*, pages 317–320, Oct 2012.

[9] K. Hoffman. *Banach spaces of analytic functions*. Prentice-Hall series in modern analysis. Prentice-Hall, 1962.

[10] J. Jugo, J. Portilla, A. Anakabe, A. Suarez, and J.M. Collantes. Closed-loop stability analysis of microwave amplifiers. *Electronics Letters*, 37(4):226–228, Feb 2001.

[11] S. Kung. A new identification and model reduction algorithm via singular value decomposition. In *Proceedings of the 12th Asilomar Conference on Circuits, Systems and Computers*, pages 705–714, 1978.

[12] Ebrahim Louarroudi. *Frequency Domain Measurement and Identification of Weakly Nonlinear Time-Periodic Systems*. PhD thesis, Vrije Universiteit Brussel (VUB), 2014.

[13] I Mallet, A. Ankabe, G. Soubercaze-Pun, and J-M Collantes. Automation of the zero-pole identification methods for the stab analysis of microwave active circuits, 2013.

[14] I. Markovsky. *Low Rank Approximation: Algorithms, Implementation, Applications*. Springer, 2012.

[15] J. P. Marmorat and M. Olivi. RARL2: a Matlab based software for $H^2$ rational approximation. http://www-sop.inria.fr/apics/RARL2/rarl2.html, 2004.

[16] Martine Olivi, Fabien Seyfert, and Jean-Paul Marmorat. Identification of microwave filters by analytic and rational $h2$ approximation. *Automatica*, 10, Februari 2013.

[17] N. Otegi, A. Anakabe, J. Pelaz, J. Collantes, and G. Soubercaze-Pun. Experimental characterization of stability margins in microwave amplifiers. *Microwave Theory and Techniques, IEEE Transactions on*, 60(12):4145–4156, Dec 2012.

[18] W. Rudin. *Real and Complex analysis*. Mc Graw-Hill, 1982.

[19] H. Sandberg, E. Mollerstedt, and Bernhardsson. Frequency-domain analysis of linear time-periodic systems. *Automatic Control, IEEE Transactions on*, 50(12):1971–1983, dec. 2005.

[20] H. Stahl. Spurious poles in padé approximation. *Journal of Computational and Applied Mathematics*, (99), 1998.

[21] A. Suarez. Check the stability: Stability analysis methods for microwave circuits. *Microwave Magazine, IEEE*, 16(5):69–90, June 2015.

[22] A. Suarez and R. Quere. *Stability analysis of nonlinear microwave circuits*. Artech House, 2002.

[23] M. van Heijningen, A. P. de Hek, F. E. van Vliet, and S. Dellier. Stability analysis and demonstration of an x-band gan power amplifier mmic. In *2016 11th European Microwave Integrated Circuits Conference (EuMIC)*, pages 221–224, Oct 2016.