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

HAL
open science

# A complementarity approach for the computation of periodic oscillations in piecewise linear systems

Valentina Sessa, Luigi Iannelli, Francesco Vasca, Vincent Acary

► To cite this version:

Valentina Sessa, Luigi Iannelli, Francesco Vasca, Vincent Acary. A complementarity approach for the computation of periodic oscillations in piecewise linear systems. Nonlinear Dynamics, 2016, 85, pp.1255-1273. 10.1007/s11071-016-2758-5. hal-02982668v2

HAL Id: hal-02982668

https://inria.hal.science/hal-02982668v2

Submitted on 29 Oct 2020

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A complementarity approach for the computation of periodic oscillations in piecewise linear systems

Valentina Sessa · Luigi Iannelli

Francesco Vasca · Vincent Acary

**Abstract** Piecewise linear (PWL) systems can exhibit quite complex behaviours. In this paper, the complementarity framework is used for computing periodic steady-state trajectories belonging to linear time-invariant systems with PWL, possibly set-valued, feedback relations. The computation of the periodic solutions is formulated in terms of a mixed quadratic complementarity problem. Suitable anchor equations are used as problem constraints in order to determine the unknown period and to fix the phase of the steady-state oscillation. The accuracy of the complementarity problem solution is shown through numerical investigations of stable and unstable oscillations exhibited by practical PWL systems: a neural oscillator, a deadzone feedback system, a stick-slip system, a repressilator and a relay feedback system.

**Keywords** Piecewise linear systems · Complementarity systems · Periodic solutions · Autonomous systems

## 1 Introduction

Piecewise linear (PWL) models can represent a wide class of practical systems and can exhibit interesting nonlinear behaviours such as periodic steady-state oscillations. In this paper, we consider the problem of the computation of periodic solutions for linear time-invariant dynamical systems whose inputs are related to the outputs through a static multi input-multi output PWL relation, possibly set-valued. That class of dynamical systems has attracted a considerable interest in the literature, and it allows to capture even complex behaviours arising in many practical systems. We can mention electrical circuits where the current-voltage characteristics of some devices can be considered as piecewise linear, such as nonlinear resistors [1], ideal diodes and switches [2,3], or mechanical systems with Coulomb friction [4]. In those applications, the steady-state behaviour is generally depicted by a periodic motion [5,6]. Since they tend to periodically oscillate also when no external excitation is applied, the period related to the periodic solution is difficult to be predicted a priori.

A considerable literature deals with the problem of computing steady-state periodic solutions and the corresponding period. The most popular time-domain

<footnote>V. Sessa (✉)
Department of Engineering, State University of Rio de Janeiro (UERJ), R. Sao Francisco Xavier, 524, Maracana, Rio de Janeiro, Brazil
e-mail: valentina.sessa@uerj.br</footnote>

<footnote>L. Iannelli · F. Vasca
Department of Engineering, University of Sannio, Piazza Roma 21, 82100 Benevento, Italy
e-mail: luigi.iannelli@unisannio.it</footnote>

<footnote>F. Vasca
e-mail: vasca@unisannio.it</footnote>

<footnote>V. Acary
INRIA, Bipop Team-Project, Inovallee de Montbonnot, 655 Avenue, Europe, 38334 Saint Ismier Cedex, France
e-mail: vincent.acary@inria.fr</footnote>
technique is the so-called shooting method [7–9]. A drawback of this method is the cost of calculating the sensitivity matrix, which is a crucial issue [10–12]. For instance, in [13, 14], in order to compute periodic stick-slip vibrations belonging to mechanical systems with discontinuous dry friction, the shooting method is applied only after the approximation of the discontinuous characteristic by a suitable smooth function. Instead, in this paper, we show that the proposed representation allows to include the discontinuity into the model. For hybrid systems that include such class, in [15] it is proposed an algorithm for the bifurcation analysis of periodic solutions. In that paper, a multi-point boundary-value problem is formulated by using the knowledge of the sequence of modes. The specific sequence is exploited also for fixing an initial guess for the periodic solution to be computed. In our paper, the periodic steady-state solution is obtained without fixing a priori the sequence of modes and without this information for choosing an initial guess.

In the frequency domain, one of the most reliable methodologies is the harmonic balance. In [16], an extension of the harmonic balance method for computing the rotary and oscillatory periodic motion of a nonlinear smooth system is proposed. Mixed time-frequency-domain approaches are presented in [17, 18]. In the first paper, a nonlinear oscillator is analysed by linearizing the system along the solution predicted by the harmonic balance technique and then by computing the Floquet's multipliers by using a time-domain numerical algorithm. Instead, in [18], the harmonic balance method is implemented together with the envelope following method in time domain. Such approach is used to compute the steady-state behaviour and the associated period of nonlinear circuits forced by two input signals with different oscillation frequencies. If one is interested in approximate results, the harmonic balance with a single harmonic leads to the describing function method [19, 20]. Such technique provides simple results about the existence of periodic solution and its parameters (the amplitude and the period), but it may fail, particularly when the system under consideration does not filter out the higher-order harmonics [21].

Complementarity models have shown to be useful for representing PWL systems [22–26] and for computing their steady-state solutions [27–29]. In [27], the authors proposed the use of a linear complementarity

approach for the evaluation of permanent oscillations in relay feedback systems with periodic exogenous signals. The autonomous case for the particular class of Lur'e systems has been analysed in [28] by representing the discretized closed-loop system in a linear complementarity form. In [29], the periodicity condition was included as a constraint of the complementarity problem, but an a priori estimation of the period was still required. Instead, in [30], the period is considered as a further unknown and it is computed together with the periodic solution by constructing a suitable mixed complementarity problem. However, such approach fails if the solution does not satisfy the assumption of being sufficiently smooth.

In this paper, we propose the use of mixed quadratic complementarity problems (MQCPs) for the computation of periodic solutions exhibited by PWL systems. The analysis here proposed extends the approach presented in [30] by allowing nonsmooth solutions and by considering a more general class of PWL systems which includes those analysed in [27–30]. In particular, the PWL system is characterized by a nonlinearity which is representable as a linear combination of saturations and relays characteristics. The effectiveness of the proposed technique is shown by investigating several examples whose steady-state solutions cannot be obtained by using the complementarity approaches previously presented in the literature.

The paper is organized as follows. In Sect. 2, we present the autonomous PWL systems of interest and we give some preliminaries about complementarity problems and the complementarity representation of a class of PWL characteristics. Section 3 shows how to formulate a mixed quadratic complementarity problem (MQCP) for computing a periodic solution together with the period. In Sect. 4, two specific cases are analysed: when the PWL system has multiple (constant and/or nonconstant) periodic solutions, and when the period of the oscillation is known. In Sect. 5, we show the effectiveness of the proposed approach by considering different steady-state behaviours in several PWL systems of practical interest. In particular, we consider a stable periodic solution in a neural oscillator, an unstable periodic solution in a deadzone feedback system, a sliding periodic solution in a stick-slip system, the periodic oscillation belonging to a repressilator and a grazing-sliding bifurcation in a relay feedback system. The paper is concluded in Sect. 6.
## 2 Preliminaries

Let us introduce first some notation and preliminary definitions. Let $\mathrm{vec}(\cdot)$ indicate the vectorization operator,<fnref n="1" /> $\mathbb{R}_+^N$ is the set of nonnegative $N$-dimensional real vectors, $I_N$ denotes the $N \times N$ identity matrix and $1_N, \infty_N$ are the $N$-dimensional vectors whose components are ones and infinity, respectively. The symbol $\otimes$ indicates the Kronecker product. The following block circulant matrix notation is adopted

$$ \Gamma_N(\mathcal{B}_1, \mathcal{B}_2) = \begin{bmatrix} \mathcal{B}_1 & 0 & \cdots & 0 & \mathcal{B}_2 \\ \mathcal{B}_2 & \mathcal{B}_1 & 0 & \cdots & 0 \\ 0 & \mathcal{B}_2 & \mathcal{B}_1 & 0 & \cdots \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & \mathcal{B}_2 & \mathcal{B}_1 \\ 0 & 0 & \cdots & 0 & \mathcal{B}_2 \end{bmatrix}, \quad (1) $$

where $\mathcal{B}_1$ and $\mathcal{B}_2$ are square matrices of the same dimension and $N$ is the number of times that $\mathcal{B}_1$ is repeated on the main diagonal.

The mixed quadratic complementarity problem (MQCP) is another important ingredient of this paper. It is defined as follows. Given a function $\varphi(z) : \mathbb{R}^r \to \mathbb{R}^r$, a square matrix $M \in \mathbb{R}^{r \times r}$, a vector $q \in \mathbb{R}^r$, and lower and upper bounds $\ell, u \in \mathbb{R}^r \cup \{-\infty, +\infty\}^r$, a MQCP is to

$$ \text{find } z \in \mathbb{R}^r, w \in \mathbb{R}_+^r, v \in \mathbb{R}_+^r \quad (2a) $$

$$ \text{s.t.} \quad \varphi(z) + Mz + q = w - v \quad (2b) $$

$$ \ell \le z \le u, \quad (z - \ell)^{\top} w = 0, \quad (u - z)^{\top} v = 0, \quad (2c) $$

where $\varphi(z)$ is a vector of quadratic forms in $z$ and the inequalities are meant componentwise. If $\ell$ and $u$ are finite, the continuity of $F(z) = \varphi(z) + Mz + q$ ensures the existence of solutions. More results about the existence of solutions for complementarity problems and variational inequality (VI) can be found in [31]. Recall

that every (nonlinear) complementarity problem as in (2) is a variational inequality (VI). When the term $\varphi(z)$ does not appear in (2a), the MQCP becomes a mixed linear complementarity problem (MLCP) [32]. Moreover, if the upper bound $u$ has all infinity components and $\ell = 0$, one gets $v = 0$ and the problem reduces to the classical linear complementarity problem with $w$ and $z$ being the usual complementarity variables [33].

The class of systems of interest can be represented by the block scheme in Fig. 1, where the dynamical system is described by

$$ \dot{x} = A_u x + B_u u + e_u \quad (3a) $$

$$ \xi = C_u x + D_u u + f_u \quad (3b) $$

$$ u = -H_u \lambda \quad (3c) $$

$$ y = H_y \xi, \quad (3d) $$

with $A_u \in \mathbb{R}^{n \times n}$, $B_u \in \mathbb{R}^{n \times p}$, $C_u \in \mathbb{R}^{p \times n}$, $D_u \in \mathbb{R}^{p \times p}$, $e_u \in \mathbb{R}^n$, $f_u \in \mathbb{R}^p$, $y = \mathrm{vec}(y_1, y_2, \dots, y_m)$, $\lambda = \mathrm{vec}(\lambda_1, \lambda_2, \dots, \lambda_m)$, $H_u \in \mathbb{R}^{m \times p}$ and $H_u \in \mathbb{R}^{p \times m}$. Each block $\mathcal{R}_i$ consists of a scalar saturation-like characteristic as in Fig. 2. Such characteristic with finite $\beta_{\mathcal{R}_i}$, $\ell_i$, $u_i \in \mathbb{R}$, $\ell_i < u_i$ and $\mu_i \ge 0$ can be represented in the following mixed linear complementarity form

$$ \mu_i \lambda_i - y_i + \beta_{\mathcal{R}_i} = w - v \quad (4a) $$

$$ \ell_i \le \lambda_i \le u_i, \quad (\lambda_i - \ell_i)w = 0, \quad (u_i - \lambda_i)v = 0, \quad (4b) $$

with $w \ge 0$ and $v \ge 0$. If $\mu_i > 0$, for any given $y_i \in \mathbb{R}$ the MLCP (4) has a unique solution $\lambda_i$ [34]. Therefore, in order to show that the MLCP (4) corresponds to the relation in Fig. 2, it is enough to verify that the feasible values of $\lambda_i$ correspond to values of $y_i$ belonging to the saturation characteristic. In particular, if $\lambda_i = \ell_i$ from (4b) it follows $v = 0$ and, since $w \ge 0$, one obtains $y_i \le \mu_i \ell_i + \beta_{\mathcal{R}_i}$, that corresponds to the lower constant saturated piece of the characteristic. If $\ell_i < \lambda_i < u_i$ from the second and third constraints in (4b) it follows $w = v = 0$, and from (4a) one obtains $y_i = \mu_i \lambda_i + \beta_{\mathcal{R}_i}$ that corresponds to the linear interval of the saturation characteristic. Finally, if $\lambda_i = u_i$, one obtains $w = 0$ and $y_i \ge \mu_i u_i + \beta_{\mathcal{R}_i}$, which corresponds to the upper saturation of the characteristic. If $\mu_i = 0$ the model (4) represents set-valued characteristics. Indeed, it is easy to verify that for $\mu_i = 0$ the model (4) represents a set-valued step function.

<footnote n="1">¹ $\mathrm{vec}\{A\}$ denotes the column vector obtained by stacking in one column all the columns of the matrix $A$. By $\mathrm{vec}(c_1, c_2, \dots, c_n)$ we will denote the column obtained by stacking all the $c_i$ columns, even when they have a different number of components.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Block diagram of a PWL system</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2** PWL saturation-like characteristic with $\mu_i > 0$</figcaption></figure>

By collecting (4) and by defining $\mathcal{R} = \text{vec}(\mathcal{R}_1, \mathcal{R}_2, \dots, \mathcal{R}_m)$, one can write the relation between $y$ and $\lambda$ as follows

$$ \text{find } \lambda \in \mathbb{R}^m, w \in \mathbb{R}_+^m, v \in \mathbb{R}_+^m \qquad (5a) $$

$$ \text{s.t. } M_{\mathcal{R}}\lambda - y + \beta_{\mathcal{R}} = w - v \qquad (5b) $$

$$ \begin{gathered} \ell_{\mathcal{R}} \le \lambda \le u_{\mathcal{R}}, (\lambda - \ell_{\mathcal{R}})^{\top} w = 0, \\ (u_{\mathcal{R}} - \lambda)^{\top} v = 0, \end{gathered} \qquad (5c) $$

with $M_{\mathcal{R}} = \text{diag}\{\mu_1, \mu_2, \dots, \mu_r\} \in \mathbb{R}^{m \times m}$, $\ell_{\mathcal{R}} = \text{vec}(\ell_1, \ell_2, \dots, \ell_m) \in \mathbb{R}^m$, $u = \text{vec}(u_1, u_2, \dots, u_m) \in \mathbb{R}^m$ and $\beta_{\mathcal{R}} = \text{vec}(\beta_{\mathcal{R}_1}, \beta_{\mathcal{R}_2}, \dots, \beta_{\mathcal{R}_m}) \in \mathbb{R}^m$.

The model (5) includes the representation of many typical PWL functions such as minimum, maximum, relay, deadzone and PWL characteristics obtained by linear combinations of these ones.

We now obtain a compact representation for the closed-loop system in Fig. 2. By substituting (3c) in (3a) and (3b), and (3b) in (3d) and by taking

$$ A = A_u \qquad (6a) $$

$$ B = B_u H_u \qquad (6b) $$

$$ C = H_y C_u \qquad (6c) $$

$$ D = H_y D_u H_u \qquad (6d) $$

$$ e = e_u \qquad (6e) $$

$$ f = H_y f_u, \qquad (6f) $$

we obtain the following dynamical system

$$ \dot{x} = Ax + B(-\lambda) + e \qquad (7a) $$

$$ y = Cx + D(-\lambda) + f \qquad (7b) $$

$$ (y, \lambda) \in \mathcal{R}, \qquad (7c) $$

with $\mathcal{R} = \text{vec}(\mathcal{R}_1, \mathcal{R}_2, \dots, \mathcal{R}_m)$ given by (5) and $(y_i, \lambda_i) \in \mathcal{R}_i, i = 1, \dots, m$.

The class of systems (7) is the one considered in this paper where $(A, B, C)$ is a minimal state space real-
ization with $A \in \mathbb{R}^{n \times n}$, $B \in \mathbb{R}^{n \times m}$, $C \in \mathbb{R}^{m \times n}$, $D \in \mathbb{R}^{m \times m}$, $e \in \mathbb{R}^n$ and $f \in \mathbb{R}^m$ all being constant and the time derivative is meant almost everywhere. It is assumed that for every initial condition $x(t_0)$ the system (7) has an absolutely continuous solution $x(t) : [t_0, +\infty) \to \mathbb{R}^n$ that satisfies (7) for almost every $t \ge t_0$. A solution $x(t)$ of the system (7) is periodic if there exists a positive real number $T$ such that $x(t+T) = x(t)$ for any $t \ge 0$. The stability concepts and definitions typically used for systems and equilibria can be also applied to periodic solutions which can be stable, asymptotically stable and unstable [35].

In the next section, it is shown how the combination of (5) with the discretized version of (7a)–(7b) allows to construct a MQCP for the computation of periodic solutions of (7).

### 3 Computation of periodic oscillations

Let $x(t)$ be a nonconstant periodic solution of (7) with unknown period, say $T$. The dynamical model (7) can be normalized with respect to the unknown period by using the time scaling $t = T\tau$, where $\tau$ is a dimensionless time variable. Then (7) can be rewritten as

$$x' = TAx + TB(-\lambda) + Te \quad (8a)$$

$$y = Cx + D(-\lambda) + f \quad (8b)$$

$$(y, \lambda) \in \mathcal{R}, \quad (8c)$$

where $x'$ is the derivative with respect to $\tau$. Any periodic solution of (7) with period $T$ will correspond to a periodic solution of (8) with period 1. Note that the right-hand side of (8a) is quadratic with respect to the unknowns $T$, $x$ and $\lambda$.

By using (5) and (8), the problem of computing the periodic solution is to

$$\begin{array}{l} \text{find } \lambda(\tau) : [0, 1] \to \mathbb{R}^m, \\ w(\tau) : [0, 1] \to \mathbb{R}_+^m, \quad v(\tau) : [0, 1] \to \mathbb{R}_+^m, \\ T \in \mathbb{R}, \quad x(\tau) : [0, 1] \to \mathbb{R}^n \end{array} \quad (9a)$$

$$\text{s. t.} \quad TB\lambda(\tau) - TAx(\tau) + x'(\tau) - Te = 0 \quad (9b)$$

$$(M_{\mathcal{R}} + D)\lambda(\tau) - Cx(\tau) - f + \beta_{\mathcal{R}} = w(\tau) - v(\tau) \quad (9c)$$

$$\begin{array}{l} \ell_{\mathcal{R}} \le \lambda(\tau) \le u_{\mathcal{R}}, \quad (\lambda(\tau) - \ell_{\mathcal{R}})^{\top} w(\tau) = 0, \\ (u_{\mathcal{R}} - \lambda(\tau))^{\top} v(\tau) = 0 \end{array} \quad (9d)$$

$$\begin{align} \lambda(0) &= \lambda(1), \quad x(0) = x(1), \notag \\ w(0) &= w(1), \quad v(0) = v(1), \tag{9e} \end{align}$$

where the constraints (9b)–(9d) must hold for any $\tau \in [0, 1]$.

The phase of a periodic solution of an autonomous system is not fixed. Then any time translation of a periodic solution provides another ‘different’ periodic solution. In other words, if the PWL system admits a periodic solution, it admits an infinite number of periodic solutions each one differing from the others by a translation in time. In order to fix the initial phase of the periodic solution, one more equation is required, which is usually called *anchor equation* [36–38]. A possible anchor equation is the phase condition proposed in [38]:

$$x_j'(\hat{\tau}) = T c_j^{\top} (Ax(\hat{\tau}) - B\lambda(\hat{\tau}) + e) = 0, \quad (10)$$

where $x_j$ is a generic $j$-th component of the state at an arbitrarily chosen time instant $\hat{\tau} \in [0, 1]$ and $c_j \in \mathbb{R}^n$ is a vector with all elements equal to zero except for the $j$-th element equal to 1. Note that the index $j$ can be chosen arbitrarily for sufficiently smooth solutions. Indeed, in case of periodic solutions $x(t)$ of class $C^1$, the time derivative of each state variable must be zero at least at one time instant $\hat{\tau} \in [0, 1]$.

In order to solve (9)–(10), one can discretize (9b) and reformulate the problem in terms of a discrete-time complementarity problem. By using the $(\theta, \gamma)$ discretization technique [39] with a sampling period $1/N$, $N$ being an integer, and by using the subscript $k$ for indicating the $k$-th sample of a variable, from (9b) one obtains

$$T\gamma B\lambda_k + T(1-\gamma)B\lambda_{k-1} - T\theta A x_k - T(1-\theta)A x_{k-1} + N x_k - N x_{k-1} - T e = 0, \quad (11)$$

where $\theta \in [0, 1]$ and $\gamma \in [0, 1]$. The constraints (9c)–(9d) can be written for each sampling time instant:

$$(M_{\mathcal{R}} + D)\lambda_k - Cx_k - f + \beta_{\mathcal{R}} = w_k - v_k \quad (12a)$$

$$\begin{align} \ell_{\mathcal{R}} &\le \lambda_k \le u_{\mathcal{R}}, \quad (\lambda_k - \ell_{\mathcal{R}})^{\top} w_k = 0, \notag \\ (u_{\mathcal{R}} - \lambda_k)^{\top} v_k &= 0, \tag{12b} \end{align}$$

with $w_k \ge 0$, $v_k \ge 0$ and $k = 1, \dots, N$. Moreover, (9e) can be rewritten as

$$\lambda_0 = \lambda_N, \quad x_0 = x_N, \quad w_0 = w_N, \quad v_0 = v_N. \quad (13)$$
Finally, one can assume that $\hat{\tau}$ in (10) is a sampling time instant. Therefore, by defining $\hat{k} = \hat{\tau}N$, the anchor equation (10) can be written as

$$T c_j^\top (Ax_{\hat{k}} - B\lambda_{\hat{k}} + e) = 0, \quad (14)$$

or, equivalently,

$$x_{\hat{k}} - x_{\hat{k}-1} = 0. \quad (15)$$

By collecting (11)-(13) for $k = 1, \dots, N$ in a matrix form we are now able to formulate the problem of computing the periodic solution (and its period) as a MQCP. In particular, the unknown vector $z$ for the MQCP is given by the samples $\lambda_k$ and $x_k$ for $k = 1, \dots, N$ and the period $T$:

$$z = \operatorname{vec}(\bar{\lambda}, \bar{x}, T) = \operatorname{vec}(\lambda_1, \dots, \lambda_N, x_1, \dots, x_N, T). \quad (16)$$

The lower and upper bounds for $T$ and for each $x_k$ with $k = 1, \dots, N$ are set as $-\infty$ and $+\infty$:

$$\ell = \operatorname{vec}(\ell_{\bar{\lambda}}, \ell_{\bar{x}}, l_T) = \operatorname{vec}(1_N \otimes \ell_{\mathcal{R}}, -\infty_{N \cdot n}, -\infty) \quad (17a)$$

$$u = \operatorname{vec}(u_{\bar{\lambda}}, u_{\bar{x}}, u_T) = \operatorname{vec}(1_N \otimes u_{\mathcal{R}}, +\infty_{N \cdot n}, +\infty). \quad (17b)$$

Given the discretization parameters $(\theta, \gamma)$ and the number of discrete samples per period $N$, by considering (16)-(17), by collecting (11)-(12) for $k = 1, \dots, N$ together with (13) and by considering (14), the MQCP can be represented in the form (2) with

$$M = \begin{bmatrix} 0 & \bar{G} & \bar{e} \\ \bar{M}_{\mathcal{R}} & \bar{C} & 0 \\ h_{\lambda} & h_{x} & 0 \end{bmatrix} \in \mathbb{R}^{(Nn+Nm+1) \times (Nn+Nm+1)} \quad (18a)$$

$$q = \operatorname{vec}(0_{Nn}, \bar{f}, \delta_q) \in \mathbb{R}^{Nn+Nm+1} \quad (18b)$$

$$\varphi(z) = \operatorname{vec}(\phi_1(z), \dots, \phi_N(z), 0_{Nm+1}), \quad (18c)$$

where $\phi_k(z) : \mathbb{R}^{Nn+Nm+1} \to \mathbb{R}^n$ are given by

$$\phi_k(z) = T\gamma B\lambda_k + T(1-\gamma)B\lambda_{k-1} \\ - T\theta Ax_k - T(1-\theta)Ax_{k-1}, \quad (19)$$

for $k = 1, \dots, N$ together with (13), and

$$\bar{G} = N\Gamma_N(I_n, -I_n) \quad (20a)$$

$$\bar{e} = -1_N \otimes e \quad (20b)$$

$$\bar{M}_{\mathcal{R}} = I_N \otimes (M_{\mathcal{R}} + D) \quad (20c)$$

$$\bar{C} = -I_N \otimes C \quad (20d)$$

$$\bar{f} = 1_N \otimes (-f + \beta_{\mathcal{R}}), \quad (20e)$$

$h_{\lambda}^{\top} \in \mathbb{R}^{Nm}$ is a vector with all zero elements, $h_{x}^{\top} \in \mathbb{R}^{Nn}$ is a vector with all zero entries but two elements $-1$ and $+1$ corresponding to the $j$-th component of $x_{\hat{k}-1}$ and $x_{\hat{k}}$, representing the anchor equation (15), and $\delta_q = 0$. The described formulation is used to compute a stable periodic oscillation in a neural oscillator and an unstable periodic solution in a deadzone feedback system in Sects. 5.1 and 5.2, respectively. Indeed, in such examples the feedback characteristic is Lipschitz continuous, then a solution of class $C^1$ is expected.

If $x(t)$ is of class $C^0$, as in the case of a relay feedback system, we have to pay attention to the fact that the time derivative can jump, so a different phase constraint should be derived. In particular, we assume that there exists a component of the system output, say $c_j^\top y$, which is equal to a value $\alpha_j$ at a time instant $\hat{\tau} \in [0, 1]$:

$$c_j^\top y = c_j^\top (Cx(\hat{\tau}) - D\lambda(\hat{\tau}) + f) = \alpha_j. \quad (21)$$

The constant $\alpha_j$ should be selected by exploiting the system structure. In particular, since we use the phase condition (21) because the trajectory $x(t)$ is only $C^0$, the discontinuity of the time derivative $\dot{x}(t)$ allows to say that there exist some $j$ and $\hat{\tau} \in [0, 1]$ such that (21) holds with $\alpha_j$ being a discontinuity point of the relation $\mathcal{R}_j$.

As in the previous case, we write (21) in discrete time as follows

$$c_j^\top (Cx_{\hat{k}} - D\lambda_{\hat{k}} + f) = \alpha_j, \quad (22)$$

with $\hat{k} = \hat{\tau}N$. The anchor equation (22) can be included in the form (18) by setting $h_{\lambda}$ and $h_{x}$ with all zero entries except for the components corresponding to $\hat{k}$ which are equal to $-c_j^\top D$ and $c_j^\top C$, respectively, and $\delta_q = c_j^\top f - \alpha_j$. In Sects. 5.3, 5.4 and 5.5, the anchor equation (22) is used to compute periodic nonsmooth oscillations in a stick-slip system, in a repressilator and in a relay feedback system.
The continuous PWL interpolation of the sequence $x_k$ obtained by solving the MQCP derived above is expected to approximate better and better the continuous-time periodic solution $x(t)$ of (7) as the number of samples per period $N$ increases. Such convergence property of the discrete-time solution is also called 'consistency' [40]. A formal proof of consistency of the solution of time-stepping methods used for solving boundary-value differential variational inequalities can be found in [41]. However, the nonsingularity condition of matrices involved in the boundary-value constraint is not valid in the case of periodic solutions.

## 4 Multiple solutions and known period

In this section, the MQCP formulated above is specified under two interesting conditions, i.e. when more periodic solutions coexist and when the period of the oscillation is known.

### 4.1 Elimination of a given periodic solution

Depending on the relation $(y, \lambda) \in \mathbb{R}$, the corresponding PWL system can present multiple periodic solutions as well as constant solutions. For instance, when $e$ and $f$ are zero and $(0, 0) \in \mathbb{R}$, the origin is also a solution (an equilibrium) of the system. Once a periodic (constant or nonconstant) solution is computed, say $\tilde{x}$, one may be interested in finding another periodic trajectory. To this aim, one can formulate an analogous MQCP by adding a further constraint which excludes the solution $\tilde{x}$ from its feasibility set. In particular, the following condition can be used:

$$ (x_{(\bullet, \hat{i})} - \tilde{x}_{(\bullet, \hat{i})})^{\top} (x_{(\bullet, \hat{i})} - \tilde{x}_{(\bullet, \hat{i})}) > 0 \quad (23) $$

for some chosen $\hat{i}$, where $x_{(\bullet, \hat{i})} \in \mathbb{R}^N$ is the vector obtained by collecting all samples of the $\hat{i}$-th component of the state vector. The condition (23) can be included in the MQCP by using the following complementarity representation:

$$ \text{find } z_{\rho} \in \mathbb{R}, w_{\rho} \in \mathbb{R}_{+}, v_{\rho} \in \mathbb{R}_{+} \quad (24a) $$

$$ \begin{array}{ll} \text{s. t.} & (x_{(\bullet, \hat{i})} - \tilde{x}_{(\bullet, \hat{i})})^{\top} (x_{(\bullet, \hat{i})} - \tilde{x}_{(\bullet, \hat{i})}) \\ & \quad - \epsilon = w_{\rho} - v_{\rho} \quad (24b) \end{array} $$

$$ \begin{array}{l} \ell_{\rho} \le z_{\rho} \le u_{\rho}, (z_{\rho} - \ell_{\rho})w_{\rho} = 0, \\ (u_{\rho} - z_{\rho})v_{\rho} = 0, \end{array} \quad (24c) $$

where $\epsilon$ is a small positive parameter, with $\ell_{\rho} = 0$ and $u_{\rho}$ equal to $+\infty$ (so that $v_{\rho} = 0$). Indeed, from the nonnegativity of the variable $w_{\rho}$, the constraints (24b) with $v_{\rho} = 0$ implies that $(x_{(\bullet, \hat{i})} - \tilde{x}_{(\bullet, \hat{i})})^{\top} (x_{(\bullet, \hat{i})} - \tilde{x}_{(\bullet, \hat{i})})$ must be strictly positive; i.e. there exists at least a sample of $x_{(\bullet, \hat{i})}$ which is different from a sample of the computed current solution $\tilde{x}_{(\bullet, \hat{i})}$.

### 4.2 MLCP for known period

A particular case of MQCP occurs when the period $T$ of the periodic solution is known. In this case, the periodic solution computation can be formulated in terms of a MLCP. Since $T$ is known, Eqs. (11) are linear with respect to the unknowns $\lambda_k$ and $x_k$. By collecting (11) and (12a) for $k = 1, \dots, N$, by defining $\ell$ and $u$ as in (17) without the last element and by using the periodicity condition (13), one obtains

$$ \begin{align} \overline{B} \overline{\lambda} + \overline{A} \overline{x} + \overline{e} &= 0 \tag{25a} \\ \overline{M}_{\mathcal{R}} \overline{\lambda} + \overline{C} \overline{x} + \overline{f} &= \overline{w} - \overline{v}, \tag{25b} \end{align} $$

with

$$ \begin{align} \overline{A} &= \Gamma_N \left( I_n - \frac{T}{N} \theta A, -I_n - \frac{T}{N} (1-\theta) A \right) \tag{26a} \\ \overline{B} &= \frac{T}{N} \Gamma_N (\gamma B, (1-\gamma)B), \tag{26b} \end{align} $$

and $\overline{e}, \overline{M}_{\mathcal{R}}, \overline{C}$ and $\overline{f}$ defined as in (20b)–(20e), respectively.

Since $T$ is known, the anchor equation (10) is not required anymore. In this case, we can define $z = \text{col}(\overline{\lambda}, \overline{x})$ and formulate the above problem as a MLCP, that is (2) without $\varphi(z)$, with the following matrices

$$ M = \begin{bmatrix} \overline{M}_{\mathcal{R}} & \overline{C} \\ \overline{B} & \overline{A} \end{bmatrix} \quad \text{and} \quad q = \begin{bmatrix} \overline{f} \\ \overline{e} \end{bmatrix}. \quad (27) $$

In the particular case of a nonsingular $\overline{A}$, one can solve (25a) for the vector $\overline{x}$ thus obtaining

$$ \overline{x} = -\overline{A}^{-1} (\overline{B} \overline{\lambda} + \overline{e}). \quad (28) $$
By substituting (28) in (25b) one obtains

$$ (\overline{M_R} - \overline{CA}^{-1}\overline{B})\overline{\lambda} + \overline{f} - \overline{CA}^{-1}\overline{e} = \overline{w} - \overline{v}, \quad (29) $$

which is in the form (2b) without $\varphi(z)$ and with $z = \overline{\lambda}$, $M = \overline{M_R} - \overline{CA}^{-1}\overline{B}$ and $q = \overline{f} - \overline{CA}^{-1}\overline{e}$. The solution $\overline{\lambda}$ of the corresponding MLCP can be used for the computation of the periodic solution $\overline{x}$ by using (28).

Note that for $T/N \to 0$, the determinant of $\overline{A}$ tends to zero. Then we can state that when a high accuracy is required, i.e. a small step size $T/N$, the formulation in (25) is preferable to (29).

In the case that the PWL system presents constant solutions, the complementarity problem can be reformulated so that the new problem has no constant solutions. Since the use of (24) makes the problem nonlinear, in order to preserve the linearity we can use an alternative approach.

As an example, let us assume that the PWL system has the origin as an equilibrium point, as discussed at the beginning of Sect. 4.1. In order to find only nonzero solutions, one can impose (without loss of generality) $x_k$ to be greater than some positive value $\varepsilon$ for some chosen $\bar{k}$. In particular, the constraint $x_{\bar{k}} \ge \varepsilon$ can be written in the complementarity form by modifying the constraints on some samples of the vector $\bar{x}$. For instance, one can set $\ell_{x_{\bar{k},1}} = \varepsilon$ and $u_{x_{\bar{k},1}} = +\infty$, where $\ell_{x_{\bar{k},1}}$, $u_{x_{\bar{k},1}}$ are the lower bound and the upper bound of the first component of $x$ at the time $\bar{k}$. Note that the approach presented above is not required in the case of an unknown period; indeed, such further constraint affects the phase condition of the solution which is already determined by the anchor equation.

## 5 Examples

In this section, the effectiveness of the proposed technique for computing period and waveform of oscillations is shown by considering PWL systems with different steady-state behaviours. We present a neural oscillator and a deadzone feedback system, which have solutions belonging to class $C^1$. The second group of examples, including a stick-slip system, a repressilator and a relay feedback system, present oscillations of class $C^0$. The computation of the periodic solution and the corresponding period is obtained by solving a MQCP with the matrices given in (18).

For all the following examples, we have computed the 'exact' steady-state solution of the discretized system by constructing the nonlinear closed-loop map. The idea for obtaining such map is to consider the solutions of the differential state equations, which model the 'modes' of the system, and to cascade them in order to link a sampled state to the next sample. It is important to highlight that the exact solution can be analytically obtained only if the sequence of modes is known, while the proposed approach is formulated without considering any hypothesis of such sequence. We numerically show that by varying the number of samples per period, i.e. by improving the resolution of the discretization, the error between the solution obtained with the MQCP and the exact one decreases. The maximum number of samples is chosen for a reasonable computational effort to perform the MATLAB code on a Intel Core i7 clocked at 2.40 GHz. The complementarity problems have been solved by using the PATH solver [32], that is a particular version of the Newton method for non-smooth problems. The solver requires an initial guess, which by default is set to zero. Since PATH provides only a local convergence, it is clear that if the initial point is close to an existing solution, one has more chances that it converges to that particular solution; e.g. one can choose as initial guess for the period the one predicted by the describing function. In the following examples, we report the initial point used in the PATH solver, whose value is chosen from a set of trivial vectors and is not related to neither to the searched solution nor to the corresponding sequence of modes. The numerical results are also compared by varying the discretization parameters.

### 5.1 Stable periodic solution in a neural oscillator

The dynamical model of a neural network consisting of $\eta$ mutually inhibiting neurons with adaptation can be represented in the following form [42]:

$$ \tau_1 \dot{\chi}_i = -\chi_i - \kappa v_i - v \sum_{j=1, j \neq i}^{\eta} u_j + e_c \quad (30a) $$

$$ \tau_2 \dot{v}_i = -v_i + u_i \quad (30b) $$

$$ u_i = \max\{0, \chi_i\}, \quad (30c) $$

with $i = 1, \dots, \eta$. Here, $\chi_i$ is the membrane potential of the i-th neuron, $v_i$ represents the degree of adaptation, $e_c$ is the total input from outside the network that
<figure><img src="image_5.png" /><figcaption>Fig. 3 Neural oscillator scheme</figcaption></figure>

is assumed positive and constant with the time, $\tau_1$ is a time constant, $\tau_2$ and $\kappa$ are the parameters that specify the time course of the adaptation, and $\nu$ indicates the strength of the inhibitory connection between the neurons. In [42], it is shown that for certain values of model parameters, the neural network has no stable equilibrium points and produces sustained oscillations. In the following, we consider a neural network with two neurons; i.e. $\eta = 2$, see Fig. 3. For neuron 1 and neuron 2, let us denote by $x_1$ and $x_3$ the state variables that represent the membrane potential, respectively, and by $x_2$ and $x_4$ the state variables for the degree of self-inhibition, respectively.

Then the neural oscillator can be represented in the form (3) with the following matrices:

$$ 
\begin{align} 
A_u &= \begin{bmatrix} -\frac{1}{\tau_1} & -\frac{\kappa}{\tau_1} & 0 & 0 \\ 0 & -\frac{1}{\tau_2} & 0 & 0 \\ 0 & 0 & -\frac{1}{\tau_1} & -\frac{\kappa}{\tau_1} \\ 0 & 0 & 0 & -\frac{1}{\tau_2} \end{bmatrix}, \nonumber \\ 
B_u &= \begin{bmatrix} 0 & -\frac{\nu}{\tau_1} \\ \frac{1}{\tau_2} & 0 \\ -\frac{\nu}{\tau_1} & 0 \\ 0 & \frac{1}{\tau_2} \end{bmatrix}, \quad e_u = \begin{bmatrix} e_c \\ 0 \\ e_c \\ 0 \end{bmatrix}, \tag{31a} \\ 
C_u &= \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix}, \nonumber \\ 
D_u &= \begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}, \quad f_u = \begin{bmatrix} 0 \\ 0 \end{bmatrix}. \tag{31b} 
\end{align} 
$$

The relation $\mathcal{R} = \text{vec}(\mathcal{R}_1, \mathcal{R}_2)$ represents the inhibitory connection between the two neurons, with $\mathcal{R}_i$ that corresponds to $\lambda_i = \max\{0, y_i\}$, $i = 1, 2$, so as indicated by (30c). Each $\mathcal{R}_i$ can be modelled as in (4) with $\mu_i = 1$, $\beta_{\mathcal{R}_i} = 0$, $\ell_i = 0$ and $u_i = +\infty$, and they can be 'linked' to the system by choosing $H_u = -I_2$ and $H_y = I_2$, in (3c) and (3d), respectively. By using (6), the neural model network is recast in form (7a) and (7b) and the feedback relation is

<figure><img src="image_6.png" /><figcaption>Fig. 4 Stable periodic steady-state oscillation of the neural oscillator projected into the $(x_1, x_3)$ plane</figcaption></figure>

represented by the complementarity form discussed above.

In [43], it is shown that by considering the parameters $\tau_1 = 0.1$ s, $\tau_2 = 0.2$ s, $\kappa = 2$, $\nu = 2$ and $e_c = 1$, the system has a locally stable periodic oscillation. Figure 4 shows the solution computed by solving the corresponding MQCP with $N = 600$, $\theta = 0.5$ and $\gamma = 0.5$ that yields a value of the period $T = 0.8973$ s. The neural system has also a constant unstable solution in [0.2, 0.2, 0.2, 0.2] that was eliminated by following the procedure presented in Sect. 4.1. The initial guess for the solver has been set as $z^0 = [x^0; T^0; \lambda^0; z^0_\rho]$, where $x^0 = [1, \dots, 1, -1, \dots, -1]$, $\lambda^0 = [1, \dots, 1, 0, \dots, 0, 1, \dots, 1]$, $T^0 = 0.897$ s, and $z^0_\rho = 0$. Let us consider the error between the exact solution and the solution computed with the complementarity approach by varying the number of samples per period and the parameters $\theta$ and $\gamma$. Since the relation $\mathcal{R}$ is a Lipschitz continuous function of $x$, the solution is expected to belong to the class $C^1$. As suggested in [3], in the case of solutions of class $C^1$, it is possible to compute the solution of the discretized system by varying both $\theta$ and $\gamma$ in order to obtain the pair which provides the best accuracy for the solution. Figure 5 shows the evolution of the maximum error in the computation of the $x_2$ trajectory, while in Fig. 6 it reported the error for the period computation.

As one could expect, the error decreases when the number of samples increases. For $\theta = 0.5$ and $\gamma$ assuming values between 0.5 and 1, the accuracy is improved, but we do not achieve order 2 when $\gamma = 0.5$. This in
<figure><img src="image_7.png" /><figcaption>Fig. 5 Neural system. Error in the steady-state solution computation obtained with different discretization techniques and with $N \in [50, 600]$: maximum error in the $x_2$ trajectory</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6 Neural system. Error in the period computation obtained with different discretization techniques and with $N \in [50, 600]$</figcaption></figure>

mainly due to the lack of regularity of the computed solution. This phenomenon is generally observed when higher-order methods are used for a solution with limited smoothness; see [4, §9.1].

## 5.2 Unstable periodic solution in a deadzone feedback system

The deadzone behaviour is typical in many real actuator systems, including mechanical connections, hydraulic

<figure><img src="image_9.png" /><figcaption>Fig. 7 Deadzone characteristic</figcaption></figure>

servo valves, piezoelectric translators and electronic circuits. A graphical representation is given in Fig. 7, and an analytical description is presented below

$$ \lambda = \begin{cases} y-1 & \text{if } y > 1 \\ 0 & \text{if } -1 \le y \le 1. \\ y+1 & \text{if } y < -1 \end{cases} \qquad (32) $$

In this section, we analyse the steady-state periodic behaviour of the dynamical system in form (3) with the following matrices:

$$ A_u = \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & -3 & -4 \end{bmatrix}, \quad B_u = \begin{bmatrix} 0 \\ 0 \\ 20 \end{bmatrix}, \quad e_u = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}, \qquad (33a) $$

$$ C_u = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}, \quad D_u = 0, \quad f_u = 0, \qquad (33b) $$

and whose feedback characteristic is the deadzone in Fig. 7. Note that $A$ is not Hurwitz, then self-induced oscillation is possible [44].

The dynamical system is a low-pass filter and the polar diagram intersects the negative real axis, so the closed-loop system is a good candidate for the application of the describing function technique [45]. The value of the approximated period is computed by considering the oscillation frequency $\omega^*$ corresponding to the point in which the polar diagram of $G(j\omega)$ (continuous line in Fig. 8) intersects the negative reciprocal of the describing function graph (dashed line in Fig. 8). We obtain $\omega^* = 1.723$ rad/s, then the corresponding period $T^* = 3.627$ s. Moreover, the orientation of the two curves predicts that the periodic oscillation is unstable.
<figure><img src="image_10.png" /><figcaption>**Fig. 8** Nyquist diagram of the transfer function of the dynamical system (continuous line) together with the representation of the negative reciprocal of the deadzone describing function (dashed line)</figcaption></figure>

The following numerical results show that the proposed complementarity approach is able to get the periodic behaviour together with the period also in the case of unstable solutions. The deadzone characteristic can be expressed as the difference of two saturation-like characteristics, $R_1$ and $R_2$. In particular, $R_1$ corresponds to $\lambda_1 = y_1$, while $R_2$ is $\lambda_2 = \text{sat}(y_2)$, that is the unitary saturation function. These two relations can be modelled as in (4) by choosing for $R_1$, $\mu_1 = 0$, $\beta_{R_1} = 0$, $\ell_1 = -\infty$ and $u_1 = +\infty$, while for $R_2$, $\mu_2 = 1$, $\beta_{R_2} = 0$, $\ell_2 = -1$ and $u_2 = +1$. Then we choose $H_u = [1, -1]$ and $H_y = \text{vec}(1, 1)$ in (3c) and (3d), respectively. By using (6), the dynamical system can be recast in form (7a) and (7b), together with the complementarity representation of feedback relation discussed above.

Figure 9 shows the periodic solution computed through the proposed technique, by solving a MQCP with matrices in (18) and, with $N = 5400$, $\theta = 0.5$ and $\gamma = 0.5$, which yields to a value of the period $T = 3.6620$ s. The initial guess for the solver has been set as $z^0 = [x^0; T^0; \lambda^0; z_\rho^0]$, where $x^0 = [1, \dots, 1, -1, \dots, -1]$, $\lambda^0 = [1, \dots, 1, -1, \dots, -1]$, $T^0 = 3.6$ s, and $z_\rho^0 = 0$. The dotted line graph corresponds to the numerical result obtained by implementing a time-stepping simulation of the initial value problem [3,46], with an initial condition sufficiently close to the periodic solution ($x_0 = [1.3, -7.848, 4]$). Of course the standard time-stepping

<figure><img src="image_11.png" /><figcaption>**Fig. 9** Unstable periodic oscillation for the deadzone feedback system: steady-state solution computed by using the MQCP approach (continuous line) and diverging time-stepping simulation (dotted line)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10** Deadzone feedback system. Error in the steady-state solution computation obtained with different discretization techniques and with $N \in [100, 5400]$: maximum error in the $x_3$ trajectory</figcaption></figure>

simulation is not able to show the unstable orbit: the proposed approach is particularly advantageous in such cases.

The error between the exact solution and the one computed with the complementarity approach has been analysed by varying the number of samples per period between 100 and 5400 and by considering values of the parameters $\theta$ and $\gamma$ in the interval $[0.5, 1]$; see Figs. 10 and 11. It is possible to prove that the trivial solution is also a solution of the system considered in this example, so the constraint (24) was added to the problem (18) to directly get the periodic solution.
<figure><img src="image_13.png" /><figcaption>Fig. 11 Deadzone feedback system. Error in the period computation obtained with different discretization techniques and with $N \in [100, 5400]$</figcaption></figure>

Qualitatively, the same comments presented for the previous example about the accuracy of the error can be repeated; see Figs. 5 and 6.

## 5.3 Sliding periodic solution in a stick-slip system

Let us consider a mass $m$ connected to a spring with elastic constant $K$, which is pulled at a constant speed $e_s$. Let $x_1$ be the elongation of the spring, and $x_2$ the mass velocity. Such system can be represented in the form (3) with

$$A_u = \begin{bmatrix} 0 & -1 \\ \frac{K}{m} & 0 \end{bmatrix}, \quad B_u = \begin{bmatrix} 0 \\ \frac{1}{m} \end{bmatrix}, \quad e_u = \begin{bmatrix} e_s \\ 0 \end{bmatrix}, \quad (34a)$$

$$C_u = \begin{bmatrix} 0 & 1 \end{bmatrix}, \quad D_u = 0, \quad f_u = 0, \quad (34b)$$

and $\mathcal{R}$ being the friction characteristic represented in Fig. 12, where $F_c$ and $F_s$ are the Coulomb friction and the stiction force, respectively. As in the previous example, such characteristic can be obtained as the combination of two saturation-like characteristics, say $\mathcal{R}_1$ and $\mathcal{R}_2$. In this case, $\mathcal{R}_1$ represents $\lambda_1 = \text{relay}(y_1)$, where relay(·) is the unitary relay characteristic, and $\mathcal{R}_2$ is $\lambda_2 = \text{sat}_{m_R}(y_2)$, where $\text{sat}_{m_R}(\cdot)$ is the unitary saturation function where the nonzero slope of the characteristic is $1/m_R$. The two relations can be modelled as in (4) by choosing for $\mathcal{R}_1, \mu_1 = 0, \beta_{\mathcal{R}_1} = 0, \ell_1 = -1$ and $u_1 = +1$, and for $\mathcal{R}_2, \mu_2 = m_R, \beta_{\mathcal{R}_2} = 0, \ell_2 = -1$ and $u_2 = +1$.

<figure><img src="image_14.png" /><figcaption>Fig. 12 Coulomb friction characteristic with Stribeck effect with null viscous friction: $y$ is the mass velocity and $\lambda$ is the friction force</figcaption></figure>

Then we choose $H_u = [F_s, -(F_s - F_c)]$ and $H_y = \text{vec}(1, 1)$ in (3c) and (3d), respectively, and we write the dynamical system in form (7) with the matrices given by (6).

Let us select $m = 1 \text{ kg}, K = 2 \text{ N/m}, F_c = 2.94 \text{ N}, F_s = 5.88 \text{ N}, m_R = 1/12$, and $e_s = 2 \text{ m/s}$. In [47], it is shown that this system can exhibit complex behaviours and that, depending on the selected integration technique, it is possible that sticking phases of the mass can be missed by the numerical simulation. Due to the discontinuity of the PWL relation at $y = 0$, a solution belonging to $C^0$ is expected. For such reason, the computation of the periodic solution and the corresponding period was obtained by solving the MQCP with matrices given in (18), using the anchor equation with the one in (22). Figure 13 shows that the complementarity approach is able to get the steady-state periodic behaviour of the system, in which it is possible to recognize the stick and the slip modes. When the trajectory is in the stick mode, the friction force increases rapidly and stops the motion. The system remains in the stick mode as long as the spring force is smaller than the stiction force. After that, the system enters in the slip mode and the mass accelerates subject to the Coulomb friction. The spring is then compressed and the motion of the mass stops again and the spring force needs to win again the stiction force. Such solu-
<figure><img src="image_15.png" /><figcaption>**Fig. 13** Periodic steady-state solution of the stick-slip system. The graphs show the elongation of the spring, $x_1$, the velocity of the mass, $x_2$, and the friction force, $\lambda$. The solution is obtained by fixing $N = 1000$, $\theta = 0.5$ and $\gamma = 1$. The corresponding value of the period is 4.9355 s</figcaption></figure>

tion is computed by solving the complementarity problem (12) with $N = 1000$, $\theta = 0.5$ and $\gamma = 1$ that yields a period $T = 4.9355$ s and by eliminating the equilibrium point in $[F_c/K, e_s]$ by following the procedure in Sect. 4.1.

The initial guess for the solver has been set as $z^0 = [x^0; T^0; \lambda^0; z^0_\rho]$, where $x^0 = [5, \dots, 5, 0, \dots, 0]$, $\lambda^0 = [(F_s - F_c), \dots, (F_s - F_c), 0, \dots, 0]$, $T^0 = 4.91$ s, and $z^0_\rho = 0$. Since the smoothness of the solution is one degree lower with respect to the previous examples, according to the discretization scheme in [4], we have fixed $\gamma = 1$ and we have varied only $\theta$. The evolution of the error shown in Figs. 14 and 15 confirms the effectiveness of the proposed approach.

## 5.4 Nonsmooth stable periodic solution in a repressilator

Let us now consider the model of a repressilator, which is a synthetic genetic regulatory network that uses a cyclic repression structure. In the literature (see, e.g., [48] and the reference therein), it is shown the existence of oscillations for a three-dimensional system in which sigmoidal regulation functions that deter-

<figure><img src="image_16.png" /><figcaption>**Fig. 14** Stick-slip system. Error in the steady-state solution computed with different discretization techniques and with $N \in [90, 1080]$; maximum error in the $x_1$ trajectory</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15** Stick-slip system. Error in the steady-state period computation obtained with different discretization techniques and with $N \in [90, 1080]$</figcaption></figure>

mine the interaction between network elements are replaced by step functions. The resulting PWL systems have been proposed as a modelling framework in biology allowing efficient simulations, as confirmed in [49] where the complementarity approach is used to compute the behaviour of a repressilator in the bacterium *Escherichia coli*. Such a repressilator consists of three genes, and it can be modelled by considering the evolution of three variables which represent concentrations of the proteins LacI, TetR and Cl. Such a repressilator model can be easily put in the form (3) with
$$
\begin{align}
A_u &= \begin{bmatrix} -a_1 & 0 & 0 \\ 0 & -a_2 & 0 \\ 0 & 0 & -a_3 \end{bmatrix}, \nonumber \\
B_u &= \begin{bmatrix} -b_1 & 0 & 0 \\ 0 & -b_2 & 0 \\ 0 & 0 & -b_3 \end{bmatrix}, \quad e_u = \begin{bmatrix} e_1 \\ e_2 \\ e_3 \end{bmatrix}, \tag{35a} \\
C_u &= \begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}, \quad D_u = 0, \quad f_u = 0, \tag{35b}
\end{align}
$$

and $\mathcal{R} = \text{vec}(\mathcal{R}_1, \mathcal{R}_2, \mathcal{R}_3)$ where

$$
\lambda_i = \begin{cases} 0 & \text{if } y_i > 1 \\ [0, 1] & \text{if } y_i = 1, \\ 1 & \text{if } y_i < 1 \end{cases} \tag{36}
$$

and $i = 1, 2, 3$. Each relation $\mathcal{R}_i$ can be modelled as the MLCP (4) with $\mu_i = 0, \beta_{\mathcal{R}_i} = -1, \ell_i = 0$ and $u_i = 1$. Then we choose $H_u = I_3$ and $H_y = -I_3$ in (3c) and (3d), respectively, and we write the dynamical system in form (7) with the matrices given by (6).

As showed in [49], the system has an unstable equi-
librium point in (1, 1, 1) and a stable nonsmooth
periodic solution. In order to compute such solu-
tion and the associated period, we solve a MQCP
with matrices in (18) replacing the anchor equation
with the one in (22). The initial guess for the solver
has been set as $z^0 = [x^0; T^0; \lambda^0; z_\rho^0]$, where $x^0 =
[3, \dots, 3, 0, \dots, 0]$, $\lambda^0 = [0, 1, 0, 0, 1, 0, \dots, 0, 1, 0]$,
$T^0 = 14.1959$ s, and $z_\rho^0 = 0$.

Figure 16 shows the time evolution of the periodic steady state for the state variables, obtained with the proposed approach. The error results in Figs. 17 and 18 show that, for this particular example, the solution is not much sensitive to the discretization parameters.

5.5 Grazing–sliding bifurcation and multiple periodic
solutions in relay feedback systems

A relay feedback system may present a great variety of
behaviours at steady state. In this section, we analyse
two examples. In the first one, the system undergoes a
grazing–sliding bifurcation, while the second example
presents multiple periodic solutions; in particular, one
solution is stable and another is unstable.

Let us consider the third-order system described by
the following matrices

<figure><img src="image_18.png" /><figcaption>**Fig. 16** Periodic steady-state solution of a repressilator network consisting of three genes. The graph shows the evolution of the variables $x_1, x_2$ and $x_3$ which represent the concentration of proteins LacI, Tetr and Cl, respectively. The solution is obtained by fixing $N = 1080$, $\theta = 0.5$ and $\gamma = 1$. The corresponding value of the period is $T = 15.1116$s</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 17** Repressilator network. Error in the steady-state solution computation obtained with different discretization techniques and with $N \in [30, 1080]$; maximum error in the $x_3$ trajectory</figcaption></figure>

$$
\begin{align}
A_u &= \begin{bmatrix} -(2\zeta\omega + \alpha) & 1 & 0 \\ -(2\zeta\omega\alpha + \omega^2) & 0 & 1 \\ -\alpha\omega^2 & 0 & 0 \end{bmatrix}, \nonumber \\
B_u &= \begin{bmatrix} \kappa \\ 2\kappa\sigma\rho \\ \kappa\rho^2 \end{bmatrix}, \quad e_u = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}, \tag{37a} \\
C_u &= \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}, \quad D_u = 0, \quad f_u = 0. \tag{37b}
\end{align}
$$
<figure><img src="image_20.png" /><figcaption>Fig. 18 Repressilator network. Error in the steady-state period computation obtained with different discretization techniques and with $N \in [30, 1080]$</figcaption></figure>

The relay characteristic can be modelled as the MLCP (4) with $M_R = 0, \beta_R = 0, \ell = -1$ and $u = 1$. We choose $H_u = 1$ and $H_y = 1$ in (3c) and (3d), respectively, and we write the dynamical system in form (7) with the matrices given by (6). In [50], it is shown that such system undergoes different type of bifurcations by varying the system parameters. We consider the grazing–sliding bifurcation that occurs when the parameter $\zeta$ varies between 0.082 and 0.098 with all the other parameters fixed at the values $\omega = 10.84, -\sigma = \alpha = \rho = \kappa = 1$. In Fig. 19, we report the periodic steady-state evolutions of the system output $x_1$ and of the corresponding relay output $\lambda$ obtained with three different values of the parameter $\zeta$. The grazing–sliding bifurcation is clearly visible. For $\zeta = 0.082$, there is a stable, symmetric trajectory with 6 pieces of sliding, which is shown in Fig. 19a. The evolution of the relay output is presented in Fig. 19b and, as expected, the value of $\lambda$ is between $-1$ and $1$ during the six time intervals corresponding to the six sliding pieces. The trajectory undergoes a grazing–sliding bifurcation and becomes a trajectory with 4 pieces of sliding for $\zeta = 0.098$ as shown in Fig. 19e, f. Figure 19c shows that for $\zeta = 0.091$ the trajectory grazes the switching surface at time instants 0.3475 and 3.243 s. Indeed, at the same time instants, the corresponding value of $\lambda$ is between $-1$ and $1$, as shown in Fig. 19d. The initial point for the solver PATH is set as $z^0 = [x^0; T^0; \lambda^0; z^0_\rho]$, where $x^0 =$

$$[1, \dots, 1, -1, \dots, -1], \lambda^0 = [-1, \dots, -1, 1, \dots, 1], T^0 = 8.642 \text{ s, and } z^0_\rho = 0.$$

Let us consider the relay feedback system with the following matrices

$$A_u = \begin{bmatrix} -2 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}, \quad B_u = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, \quad e_u = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}, \qquad (38a)$$

$$C_u = \begin{bmatrix} 0.02 & 0.3 & 1 \end{bmatrix}, \quad D_u = 0, \quad f_u = 0. \qquad (38b)$$

As in Sect. 5.2, let us use the describing function technique for predicting the existence of periodic solutions. Figure 20 shows the Nyquist diagram of system (38) together with the negative reciprocal of the relay describing function, which corresponds to the negative real axis. The Nyquist diagram intersects the negative real axis in two points in which $\omega$ assumes the values 1.65 and 4 rad/s. This means that the closed-loop system has two periodic solutions with different period. In particular, the system has an unstable periodic solution of period $T = 1.570$ s and a stable periodic solution of period $T = 3.808$ s.

By using the complementarity representation for the relay characteristic given above, we can write the dynamical system in form (7) with the matrices given by (6).

The unstable solution in Fig. 21 and the corresponding period are computed by solving the MQCP in (18) and by setting $N = 400, \gamma = 1$ and $\theta = 1$. The computed value of the period is 1.8030 s. The initial point for the solver PATH is set as $z^0 = [x^0; T^0; \lambda^0]$, where $x^0 = [0, \dots, 0, 1, \dots, 1]$, $\lambda^0 = [1, \dots, 1, 0, \dots, 0]$, and $T^0 = 2$ s. For the dynamical system under consideration, the constant trivial solution is also a solution. By using the arguments in Sect. 4, suitable constraints were added to the MQCP in order to eliminate both the computed periodic unstable solution and the trivial one. Then the periodic stable steady-state solution in Fig. 22 was obtained, whose period is 3.5645 s. We use the same initial point, except that we added the initial value for the variables $z_{\rho_1}$ and $z_{\rho_2}$, which was set to zero.

## 6 Conclusions

The mixed quadratic complementarity approach has been proposed for the computation of periodic solu-
<figure><img src="image_21.png" /></figure>

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

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

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

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

<figure><img src="image_26.png" /><figcaption>Fig. 19 Relay feedback system: effect of a grazing–sliding bifurcation. Periodic steady-state trajectory of the variable $x_1$: **a** for $\zeta = 0.082$ it shows six pieces of sliding, **c** and for $\zeta = 0.091$ the trajectory grazes the switching surface at time instants 0.3475</figcaption></figure>

and 3.243 s, **e** and for $\zeta = 0.098$ the trajectory has four sections of sliding. Periodic steady-state evolution of the relay output: **b** for $\zeta = 0.082$, **d** for $\zeta = 0.091$ and **f** for $\zeta = 0.098$
<figure><img src="image_27.png" /><figcaption>**Fig. 20** Nyquist diagram of the transfer function of the dynamical system (38) (*continuous line*) together with the representation of the negative reciprocal of the relay describing function (*dashed line*)</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 21** Unstable periodic steady-state solution of the relay feedback system in (38). The solution is obtained by fixing $N = 400$, $\theta = 1$ and $\gamma = 1$. The corresponding value of the period is 1.8030 s</figcaption></figure>

tions in a class of piecewise linear (PWL) systems, which consists of linear time-invariant systems with PWL feedback relations representable as linear combinations of saturation-like or step characteristics. For many practical systems belonging to this class, the period and the shape of the oscillation are difficult to be predicted, then phase conditions acting as anchor equations for the periodic solution have been added to the complementarity problem. It has been showed how to build a mixed quadratic complementarity problem (MQCP) which models the discretized closed-loop system together with the periodicity constraint and the

<figure><img src="image_29.png" /><figcaption>**Fig. 22** Stable periodic steady-state solution of the relay feedback system in (38). The solution is obtained by fixing $N = 400$, $\theta = 1$ and $\gamma = 1$. The corresponding value of the period is 3.5645 s</figcaption></figure>

phase condition. The solution of such MQCP gives, as outcome, the periodic solution and the corresponding period. The effectiveness of the proposed approach has been tested by considering several practical PWL systems with different steady-state periodic behaviours: a stable periodic solution in a neural oscillator, an unstable oscillation for a deadzone feedback system, a sliding orbit in a stick-slip system, a periodic solution belonging to a repressilator, a grazing-sliding bifurcation in a relay feedback system and coexisting periodic solutions in a relay feedback system. The simulations, implemented by using different discretization schemes and sampling frequencies, showed the good accuracy of the proposed approach with respect to the analytical solution, computed by assuming the a priori knowledge of the modes sequence. Future work will investigate the use of the proposed tool for the computation of periodic solutions in a continuation framework.

## References

1. Stern, T.E.: Piecewise-linear network theory. Technical report, PhD thesis, MIT, Research Laboratory of Electronics, Cambridge, MA, USA (1956)

2. Camlibel, M.K., Heemels, W.P.M.H., van der Schaft, A.J., Schumacher, J.M.: Switched networks and complementarity. IEEE Trans. Circuits Syst. I Reg. Pap. 50(8), 1036–1046 (2003)

3. Acary, V., Bonnefon, O., Brogliato, B.: Time-stepping numerical simulation of switched circuits within the non-smooth dynamical systems approach. IEEE Trans. Comput. Aided Des. Integr. Circuits Syst. 29(7), 1042–1055 (2010)
4. Acary, V., Brogliato, B.: Numerical Methods for Nonsmooth Dynamical Systems. Springer, Berlin (2008)

5. Melin, J., Hultgren, A., Lindstrom, T.: Two types of limit cycles of a resonant converter modelled by three-dimensional system. Nonlinear Anal. Hybrid Syst. 2(4), 1275-1286 (2008)

6. Yang, X.-S., Chen, G.: Limit cycle and chaotic invariant sets in autonomous hybrid planar systems. Nonlinear Anal. Hybrid Syst. 2(3), 952-957 (2008)

7. Aprille, T., Trick, T.: A computer algorithm to determine the steady-state response of nonlinear oscillators. IEEE Trans. Circuit Theory 19(4), 354-360 (1972)

8. Flieller, D., Riedinger, P., Louis, J.P.: Computation and stability of limit cycles in hybrid systems. Nonlinear Anal. Theory Methods Appl. 64(2), 352-367 (2006)

9. Li, D., Xu, J.: A new method to determine the periodic orbit of a nonlinear dynamic systems and its period. Eng. Comput. 20(4), 316-322 (2005)

10. Leine, R.I., Nijmeijer, H.: Dynamics and Bifurcations of Non-Smooth Mechanical Systems. Springer, London (2004)

11. Theodosiou, C., Pournaras, A., Natsiavas, S.: On periodic steady state response and stability of Filippov-type mechanical models. Nonlinear Dyn. 66(3), 355-376 (2011)

12. Bizzarri, F., Brambilla, A., Gajani, G.S.: Steady state computation and noise analysis of analog mixed signal circuits. IEEE Trans. Circuits Syst. I Reg. Pap. 59(3), 541-554 (2012)

13. Leine, R.I., Van Campen, D.H., De Kraker, A.: Stick-slip vibrations induced by alternate friction models. Nonlinear Dyn. 16(1), 41-54 (1998)

14. Van De Vrande, B.L., Van Campen, D.H., De Kraker, A.: An approximate analysis of dry-friction-induced stick-slip vibrations by smoothing procedure. Nonlinear Dyn. 19(2), 157-169 (1999)

15. Thota, P., Dankowicz, H.: TC-HAT (TC): a novel toolbox for the continuation of periodic trajectories in hybrid dynamical systems. SIAM J. Appl. Dyn. Syst. 7(4), 1283-1322 (2008)

16. Lu, C.-J., Lin, Y.-M.: A modified incremental harmonic balance method for rotary periodic motions. Nonlinear Dyn. 66(4), 781-788 (2011)

17. Bonnin, M., Gilli, M., Civalleri, P.P.: A mixed time-frequency-domain approach for the analysis of a hysteretic oscillator. IEEE Trans. Circuits Syst. II Exp. Br. 52(9), 525-529 (2005)

18. Brambilla, A., Gruosso, G., Storti Gajani, G.: MTFS mixed time-frequency method for the steady-state analysis of almost-periodic nonlinear circuits. IEEE Trans. Comput. Aided Des. Integr. Circuits Syst. 31(9), 1346-1355 (2012)

19. Duarte, F.B., Tenreiro Machado, J.: Fractional describing function of systems with Coulomb friction. Nonlinear Dyn. 56(4), 381-387 (2009)

20. Huang, Y.-J., Wang, Y.-J.: Steady-state analysis for a class of sliding mode controlled systems using describing function method. Nonlinear Dyn. 30(3), 223-241 (2002)

21. Engelberg, S.: Limitations of the describing function for limit cycle prediction. IEEE Trans. Autom. Control 47(11), 1887-1890 (2002)

22. Vasca, F., Iannelli, L., Camlibel, M.K.: A new perspective for modeling power electronics converters: complementarity framework. IEEE Trans. Power Electron. 24(2), 456-468 (2009)

23. Schumacher, J.M.: Complementarity systems in optimization. Math. Program. 101(1), 263-296 (2004)

24. Sessa, V., Iannelli, L., Vasca, F.: A complementarity model for closed-loop power converters. IEEE Trans. Power Electron, 29(12), 6821-6835 (2014)

25. van der Schaft, A.J., Schumacher, J.M.: Complementarity modeling of hybrid systems. IEEE Trans. Autom. Control 43, 483-490 (1998)

26. Heemels, W.P.M.H., De Schutter, B., Bemporad, A.: Equivalence of hybrid dynamical model. Automatica 37(7), 1085-1091 (2001)

27. Iannelli, L., Vasca, F.: Computation of limit cycles and forced oscillations in discrete-time piecewise linear feedback systems through a complementarity approach. In: 47th IEEE Conference on Decision and Control, pp. 1169-1174, Cancun, Mexico (2008)

28. Iannelli, L., Vasca, F., Sessa, V.: Computation of limit cycles in Lur'e systems. In: American Control Conference, pp. 1402-1407, San Francisco, CA, USA (2011)

29. Sessa, V., Iannelli, L., Vasca, F.: Mixed linear complementarity problems for the analysis of limit cycles in piecewise linear systems. In: The 51st IEEE Conference on Decision and Control, pp. 1023-1028, Maui, Hawaii, USA (2012)

30. Sessa, V., Iannelli, L., Acary, V., Brogliato, B., Vasca, F.: Computing period and shape of oscillations in piecewise linear Lur'e systems: a complementarity approach. In: The 52nd IEEE Conference on Decision and Control, pp. 4680-4685, Florence, Italy (2013)

31. Facchinei, F., Pang, J.S.: Finite-Dimensional Variational Inequalities and Complementarity Problems. Operations Research. Springer, New York (2003)

32. Dirkse, S.P., Ferris, M.C.: The PATH solver: a non-monotone stabilization scheme for mixed complementarity problems. Optim. Methods Softw. 5, 123-156 (1995)

33. Cottle, R., Pang, J., Stone, R.: The Linear Complementarity Problem, 2nd edn. Academic Press, Cambridge (2009)

34. Gowda, M.S., Pang, J.: Stability analysis of variational inequalities and nonlinear complementarity problems, via the mixed linear complementarity problem and degree theory. Math. Oper. Res. 19(4), 831-879 (1994)

35. Farkas, M.: Periodic Motions. Series in Applied Mathematical Sciences edn. Springer, New York (1994)

36. Doedel, E.J.: AUTO: a program for the automatic bifurcation analysis of autonomous systems. Congr. Numer. 30, 1265-1284 (1981)

37. Nayfeh, A.H., Balachandran, B.: Applied Nonlinear Dynamics: Analytical, Computational and Experimental Methods. Wiley, Chichester (1995)

38. Seydel, R.U.: Practical Bifurcation and Stability Analysis. Interdisciplinary Applied Mathematics, 3rd edn. Springer, New York (1988)

39. Acary, V.: Time-stepping via complementarity. In: Vasca, F., Iannelli, L. (eds.) Dynamics and Control of Switched Electroinic Systems, Advances in Industrial Control. Springer, London (2012)

40. Han, L., Tiwari, A., Camlibel, M.K., Pang, J.-S.: Convergence of time-stepping for passive and extended linear complementarity systems. SIAM J. Numer. Anal. 47(5), 3768-3796 (2009)

41. Pang, J.-S., Stewart, D.E.: Differential variational inequalities. Math. Program. 113, 345-424 (2008)
42. Matsuoka, K.: Sustained oscillations generated by mutually inhibiting neurons with adaptation. *Biol. Cybern.* **52**, 367-376 (1985)

43. Goncalves, J.M.: Region of stability for limit cycles in piecewise linear systems. *IEEE Trans. Autom. Control* **50**(11), 1877–1882 (2005)

44. Hu, T., Thibodeau, T., Teel, T.R.: A unified Lyapunov approach to analysis of oscillations and stability for systems with piecewise linear elements. *IEEE Trans. Autom. Control* **55**(12), 2864–2869 (2010)

45. Khalil, H.K.: *Nonlinear Systems*. Prentice Hall, New Jersey (2002)

46. Heemels, W.P.M.H., Camlibel, M.K., Schumacher, J.M.: A time-stepping method for relay systems. In: 39th IEEE Conference on Decision and Control, pp. 4461–4466, Sydney, Australia (2000)

47. Astrom, K.J., Canudas De Wit, C.: Revisiting the LuGre friction model. *IEEE Control Syst. Mag.* **28**(6), 101–114 (2008)

48. Tonnelier, A.: Cyclic negative feedback systems: what is the chance of oscillation? *Bull. Math. Biol.* **76**(5), 1155–1193 (2014)

49. Acary, V., de Jong, H., Brogliato, B.: Numerical simulation of piecewise-linear models of gene regulatory networks using complementarity systems. *Phys. D Nonlinear Phenom.* **269**, 103–119 (2014)

50. Bernardo, M., Budd, C., Champneys, A.R., Kowalczyk, P.: Piecewise-Smooth Dynamical Systems, vol. 163. Springer, London (2008)