Bifurcation analysis for non-local design of a hybrid observer for the impulsive Goodwin’s oscillator

The impulsive Goodwin’s oscillator is a mathematical model capturing the dynamics arising in a closed-loop system, where a third-order linear time-invariant plant is subject to an intrinsic pulse-modulated feedback. Originally, the model was motivated by pulsatile regulation in endocrine systems but also has other potential applications. The asymptotic estimation of the hybrid state of the impulsive Goodwin’s oscillator is considered in this paper. A hybrid observer makes use of the continuous plant output to correct the estimates of the state vector through two output error feedbacks: a continuous and a discrete one. When the hybrid state estimation error is zero, the observer is in a synchronous mode characterized by the firings of the impulses in the observer feedback and those of the plant occurring simultaneously. The synchronous mode thus corresponds to an equilibrium point of the hybrid state error dynamics. To guarantee asymptotic convergence of the observer to the synchronous mode, the basin of attraction of the equilibrium has to include all feasible initial deviations of the state estimates. To guarantee the above properties, a numerical observer design approach based on bifurcation analysis of a discrete map capturing the observer state transitions from one impulse firing to another is proposed and its efficacy is demonstrated in simulation.


Introduction
The impulsive Goodwin's oscillator is a mathematical model comprised of a linear third-order time-invariant continuous plant controlled by a pulse-modulated feedback, [1,2]. The closed-loop dynamics are hybrid and lack equilibria due to the form of the amplitude modulation function. On the one hand, the model generalizes the classical continuous Goodwin's oscillator [3][4][5][6][7][8][9] and provides a modeling framework to describe numerous oscillatory biological control systems exploiting neurally implemented impulsive feedback. Exactly as its predecessor, the impulsive Goodwin's oscillator possesses positivity of the state variables. On the other hand, it presents an analytically tractable example of a hybrid oscillator with rich nonlinear dynamics including, e.g., high-periodic solutions and determin-istic chaos. At the same time, and in contrast with the continuous counterpart, the considered impulsive oscillator is validated through system identification on human endocrine data with respect to testosterone regulation in the male, [10]. The latter fact secures the utility of the model. A structurally similar to the impulsive Goodwin's oscillator model has been more recently proposed to describe the pharmacokinetics of levodopa, the mainstay drug in Parkinson's disease, subject to intermittently interrupted gastric emptying.
The impulsive feedback signal of the Goodwin's oscillator is intrinsic to the closed loop and essentially cannot be measured in biomedical applications, which property poses an interesting and practically meaningful problem for estimating the full hybrid state of a system from only continuous measurements. In endocrinology, the problem of reconstructing the episodic firings of a pulsatile feedback is standard and usually handled by applying deconvolution algorithms under assumptions on the signal shape of the pulse; see, e.g., [11,12]. While deconvolution algorithms are widely used in practice due to readily available software, being applied to pulsatile endocrine regulation, they basically ignore the pulse-modulated feedback nature of the underlying system and are not suitable for online implementation. The deconvolution per se (i.e., blind deconvolution) is not a well-posed problem and needs regularization. Further, all the deconvolution methods used in practice enforce, in some way, solution sparsity and merge the pulses near to each other to one. Faced with a signal shape that deviates much from linear dynamics, the deconvolution methods allow for multiple pulses within a short time interval. This results in a deceptively insignificant output modeling error but yet does not properly capture the impulsive feedback characteristics. Therefore, it is motivated to consider output error feedback observer algorithms for the reconstruction of the impulsive feedback action along with the continuous states from a continuous output since they project measured data onto an analytically defined mathematical model.
Notice that in most of the studies covering state reconstruction in hybrid system, the discrete states (or switching events) are assumed to be known, as they typically belong to the engineered part of the system, as, e.g., in automotive applications [13]. Because control and observation over computer network become ubiquitous, the opposite problem of estimating the states of a continuous system from irregular discrete measure-ments is considered much more broadly and actively, e.g., [14].
Observers for hybrid systems that utilize only continuous measurements are rarely considered in the literature and limited to specific plant types. Since the impulsive Goodwin's oscillator only exhibits periodic and non-periodic (i.e., quasiperiodic and chaotic) oscillations, the observer design problem can be formulated as synchronization of the plant states with those of the observer. This approach has been proposed in [15] and does not formally demand observability of the hybrid plant solution. In a certain way, the idea is similar to that of the observers for unobservable but detectable linear time-invariant systems, where the state estimates that cannot be corrected by output error feedback converge due to the stable intrinsic dynamics. Similarly, the dynamics governing the observed periodical solution to the impulsive Goodwin's oscillator have to be orbitally stable for the observer to be functional. Further, since the oscillator dynamics are highly nonlinear, the state estimation problem is considered with respect to particular solutions of the observed system, whose characteristics are assumed to be known, but not the initial conditions. This can be compared to the situation with the Luenberger observer, where the model of the plant is available and the state estimation error due to the mismatch in the initial conditions on the plant and the observer asymptotically vanishes.
The observer design problem for the impulsive Goodwin's oscillator consists of assigning desired properties to a discrete map that captures the observer state transitions from one impulse firing to another through manipulating the degrees of freedom provided by the observer structure. A simplest observer structure includes just one output error feedback matrix for the continuous state estimates, as proposed in [15]. Then, it is proven that there is a sufficiently high observer gain that guarantees asymptotic convergence of all the components of the hybrid state estimate.
Since the simplest observer structure suffers from slow convergence in the discrete state estimate (i.e., the impulsive feedback firing time), an additional observer gain from the output error to the discrete state estimate is introduced in [16]. In fact, even with a slightly unstable (physiologically infeasible) linear part of the model, the hybrid state estimation error can be asymptotically stabilized solely by the gain to the discrete state estimate; see [17]. With two nonzero gains, the observer design becomes an intricate highly nonlin-ear problem of balancing the gain to the continuous state variable estimates against that to the discrete state estimate in order to achieve a desirable overall performance. With high continuous gain, the estimation error in the continuous states converges fast which implies that the (continuous) output estimation error is virtually zero and does not provide correction to the discrete state. When the discrete observer gain is set too high, the hybrid state estimate convergence is disrupted by so-called overjumps that result in postponing the synchronization of the plant and observer by one or, more often, several periods of the observed plant solution.
Control and observer design for nonlinear and hybrid systems is heavily focused on analytical techniques and global stability properties of the designed structure. This is in a sharp contrast with the engineering practices that are based on simulation of complex plants and tuning of control and estimation algorithms through more or less formal optimization algorithms. Bifurcation analysis is an established tool for discerning and studying in depth the mechanisms of complex dynamics of a given system. Applying bifurcation analysis to design of nonlinear control is not a novel idea but more accepted and utilized in certain specific areas of engineering that are challenged with nonlinear dynamical phenomena. Flight control [18], chemical reactor control [19], and closed-loop control of drug delivery [20] are just few examples of nonlinear control applications, where bifurcation analysis and continuation techniques serve design purposes well. However, analysis of nonlinear dynamics becomes computationally heavy in higher dimensions which fact limits the number of controller/observer design degrees of freedom. Further, the structure of a control or estimation algorithm to be analyzed has to be specified beforehand.
The problem considered in this paper is how to select the design parameters of a hybrid state observer that render appropriate state estimation performance. Since the observer is based on a model describing the plant, the achievable state reconstruction performance is restricted by the complex nonlinear plant dynamics. Bifurcation analysis is used to discover and observe these limitations in the proposed observer design algorithm.
Novelty of the proposed approach is ensured by the fact that bifurcation analysis for hybrid observer design has not been covered before in the literature. The main contributions of this study are three and stand as follows: -Firstly, to enable a bifurcation analysis, an approach to determine the periodic motions and investigating their local stability is proposed. The problem of periodic motion evaluation is reduced to solving a system of transcendental equations with respect to the fixed points of the periodic orbit of the discrete mapping. This allows for the evaluation of stable as well as unstable cycles. -Secondly, a complete bifurcation analysis of a hybrid observer for the impulsive Goodwin's oscillator is performed elucidating the complex dynamics phenomena that arise in it. -Finally, a systematic design procedure that guaranties the observer function for all feasible initial conditions of the hybrid state estimate is proposed.
The rest of the paper is structured as follows. First, the equations of the impulsive Goodwin's oscillator are recapitulated and the main dynamical properties of it are summarized. To reconstruct the hybrid state vector from the continuous output, an observer structure equipped with two output error feedback gains, one correcting the continuous state variables and another adjusting the discrete one, is introduced. Further, a discrete map describing the hybrid dynamics of the observer is derived and studied by bifurcation analysis highlighting the distinctive nonlinear phenomena. Based on those findings, a systematic observer design algorithm is proposed and illustrated by numerical examples.

The impulsive Goodwin's oscillator
The impulsive version of Goodwin's oscillator has been introduced in [1] to portray the hybrid dynamics of testosterone (Te) regulation arising due to the pulse-modulated feedback implemented by the hypothalamus. Pulses of gonadotropin releasing hormone (GnRH) are secreted from the hypothalamus to stimulate the release of luteinizing hormone (LH) from the pituitary gland that, in its turn, stimulates the release of testosterone (Te) from the testes. Then, Te inhibits the secretion of GnRH through a negative feedback by lowering the amplitude and frequency of the pulses.
The continuous part of the impulsive Goodwin's oscillator is given bẏ where x 1 (t) undergoes jumps at the time instants t n , n > 0: In terms of Te regulation in the male, x 1 is interpreted as the concentration of GnRH, x 2 as the concentration of LH, and x 3 as the concentration of Te. The positive parameters b 1 , b 2 and b 3 denote the reciprocal half-life times of the hormones, and g 1 and g 2 are the secretion rate factors of LH and Te, respectively. The nonlinear functions F(·) and Φ(·) are positive and bounded; F(·) is decreasing and Φ(·) is increasing. In particular, the positivity of F(·) explains why there are no equilibria in the impulsive Goodwin's oscillator. The superscripts "− and "+ denote the left-and rightside limits, respectively. The concentration of GnRH is subject to jumps due to the pulse-modulated frequency and amplitude feedback from Te implemented via the hypothalamic neurons. The positivity of Φ(·) corresponds to the necessity for the neuron to have a refractory period between firings.
The impulsive hormonal feedback modulation functions F(·) and Φ(·) are assumed to be Hill functions of the form: where Φ 1 , Φ 2 , F 1 , F 2 , h are positive constants and β is the order of the Hill functions. For biological feasibility of the model, the following bounds apply: With the state vector x = [x 1 x 2 x 3 ] , model Eq. (1) can be rewritten aṡ Thus, y(t) comprises the measurable serum concentrations of LH and Te, and z(t) is the concentration of Te that drives the impulsive feedback. Since the equalities C B = 0, L B = 0 hold, the outputs y(t) and z(t) are continuous. Notice that all eigenvalues of the matrix A are real and negative. This property corresponds to the fact that all organic molecules eventually decay. In addition, the matrix A is also Metzler. This implies that the positive octant R 3 + is invariant with respect to Eq. (3) and the continuous part of Eq. (3) is a positive system. System positivity is an important property in mathematical modeling of biological systems, in particular those whose state variables stand for concentrations of chemicals. It can be checked that the pair (A, L) is observable and the fact immediately follows from the cascaded structure of Eq. (1).
As shown in [1], any solution of model Eq. (1) is bounded and the considered impulsive feedback system does not possess an equilibrium point. For biologically motivated parameter values, dynamical system Eq. (3) exhibits a stable periodic solution with one or two pulses in the period, i.e., a stable 1-cycle, or a stable 2cycle. However, as demonstrated in [21], more complicated phenomena are also possible, including perioddoubling transition to chaos, period-adding phenomena and quasiperiodicity. Importantly, the impulsive Goodwin's oscillator is currently the only endocrine feedback regulation model that is validated on human Te and LH data [10].
The pulses of amplitude λ n produced by the impulsive feedback at the firing times t n are immeasurable thus motivating their model-based reconstruction. Therefore, the state observation problem at hand consists of producing the estimates (t n ,λ n ) of the firing times t n and the amplitudes λ n of the pulses. In the endocrine closed-loop system, the GnRH pulses cannot be measured for ethical reasons as they are secreted deep in the brain and the hormone has a short halflife time. Note that, with a completely known impulse sequence, estimates of the continuous state vector x(t) can be obtained by conventional state estimation techniques.

Observer structure
In order to estimate the state vector of system Eq. (3), a hybrid observer is introduced in [16] aṡ witĥ where K ∈ R 3×2 is the continuous feedback matrix and k d is the gain in the discrete part of the observer. The gains K , k d represent the design degrees of freedom. The discrete dynamics of Eq. (4) are directly impacted by the scalar gain k d , while K controls the convergence of the continuous estimatex(t) in between the firings of the pulse-modulated feedback.
Notably, even with zero value of either gain, the observer is still able to converge to the hybrid state vector of the plant. With k d = 0, the observer is known to suffer from slow convergence; see [15]. Since the pair (A, L) is observable, the gain K assigns the convergence rate of the continuous state estimation error. As the matrix A is Hurwitz stable, the continuous gain is not necessary. In fact, even the gain matrices K that result in (slightly) non-Hurwitz matrices can still produce stable hybrid dynamics through the stabilizing action of k d , as the continuous and discrete observer dynamics interact.

Synchronous mode
To completely characterize the hybrid system dynamics of Eq. (3), not only the continuous state vector x(t), but also the sequence of time instants t n should be taken into account.
Let (x(t), t n ) be a solution of Eq. (3) that undergoes jumps at time instants t n , n 0: and denote x n = x(t − n ). Suppose that the hybrid state (x(t), t n ) is reconstructed by means of observer Eqs. (4)- (6). Without loss of generality, assume that t 0 t 0 < t 1 . Obviously, the solution (x(t),t n ) of observer equations (4)- (6) subject to the initial conditionst 0 = t 0 , is called a synchronous mode of observer Eqs. (4)- (6) with respect to the solution (x(t), t n ) of plant Eq. (3); see [15]. Thus, a synchronous mode corresponds to the zero dynamics of the hybrid observer.
Since the continuous plant output y(t) is the only source of information about the hybrid plant state, the result below is instrumental.
The following conditions are equivalent Conditions in Proposition 1 ensure that reducing the plant output estimation error to zero brings the observer into a synchronous mode and, therefore, solves the hybrid state estimation problem.
To ensure practical usefulness of the observer, stability properties of the synchronous mode have to be investigated. Indeed, once perturbed, the observer in a synchronous mode has to rapidly return into it. A synchronous mode with respect to (x(t), t n ) is called locally asymptotically stable (see [15]) if, for any solu- Observing a periodic solution (x(t), t n ) of system Eq. (3) with exactly m 1 pulses within the least period of the solution is of special interest. Such a solution is called m-cycle. Then, x n+m ≡ x n , λ n+m ≡ λ n , T n+m ≡ T n . Logically, the synchronous mode (x(t),t n ) is also an m-cycle withx n+m ≡x n ≡ x n , λ n+m ≡λ n ≡ λ n ,T n+m ≡T k ≡ T n . Denote the mcycles of the plant and the observer by C m (x(t), t n ) and C m (x(t),t n ) , respectively.
Consider the set It was proved in [1] that, for a trajectory starting in Z (i.e., x(0) ∈ Z ), the following inequalities hold , consider also the set Obviously, T Σ is the least period of the m-cycle being estimated. Therefore, the space D = Z × T (where × is a Cartesian product) constitutes the (phase) subspace of feasible initial conditions for the observer.
Introduce the continuous state estimation error For all values of t, where r (t) is continuous, it satisfies the differential equatioṅ The function r (t) undergoes jumps at the points t =t n . Thus, when evolving in t, r (t) is subject to jumps induced both by the impulsive feedback in the plant and that of the observer. In general, the continuous state estimation error r (t) can belong to: -an equilibrium point corresponding to a synchronous mode, -a periodic attractor that contains solutions with asynchronous firing of impulses in the plant and the observer, -a quasiperiodic attractor, or -a chaotic attractor.
Of the alternatives above, only the first one is acceptable when an accurate hybrid state estimate is desired. The second option in the list is an undesirable behavior featuring, for a certain values of K and k d , an p-periodic asynchronous mode withx n+ p ≡x n ,t n = t n , where p 1 is some integer (Fig. 1). If such asynchronous mode is locally stable, then, for certain initial conditions, the observer might enter this mode, instead of a synchronous one.
To include the whole space of feasible initial conditions D into the domain of attraction of the synchronous mode, the gain coefficients K and k d have to be selected so that there are no other periodic modes of the observer, distinct from the synchronous one, or these modes are unstable. The conditions for the existence and stability of a periodic asynchronous mode are provided in the next section. Further, by choice of K and k d , the synchronous mode has to be rendered locally asymptotically stable and attractive on D with a suitable convergence rate. This constitutes the observer design problem at hand.

Discrete mapping
Recall from [17] the discrete mapping capturing the propagation of the continuous observer states through the discrete cumulative (plant and observer) sequence of the feedback firing instants Q : Hence, each point of the observer hybrid state (x n ,t n ) belongs to one of the sets S k,s , i.e., to each (x n ,t n ) one can uniquely match two points (x k , t k ) and (x s , t s ) of the observed system (if k = s, these points coincide) such that t k t n < t k+1 , t s t n + Φ(α(x n ,t n )) < t s+1 . Now mapping Eq. (10) can be written as As proven in [17], the mapping Q(ζ, θ ) is smooth provided the functions Φ(·) and F(·) have continuous derivatives. Thus, the discrete observer dynamics admit linearization in the vicinity of (x n ,t n ), where local stability properties of observer solutions can be analyzed.
Disregarding for a moment the task of state observation, one can look upon observer Eqs. (4)-(6) as a hybrid dynamical system driven by the plant output y(t). Then, for a given plant model, each solution of the observer equations is completely defined by the initial conditions on the observerx(t − 0 ), the observer gains K , k d , and the time evolution of the plant output y(t). Therefore, in view of synchronous and asynchronous observer modes (e.g., as in Fig. 1), it is now possible to obtain general conditions for existence of a periodical solution in the observer.
Consider an m-cycle C m (x(t), t n ) of plant Eq. (3) that is characterized by a sequence of the jump points x n and the intervals T n = t n+1 − t n , n = 0, . . . , m − 1, T 0 = T m between the subsequent firing times. Assume that there exists a p-cycle C p (x(t),t n ) of observer Eqs. (4)-(6), such that each point (x n ,t n ) (n = 0, . . . , p−1) corresponds to two points (x k n , t k n ) and (x s n , t s n ) of the plant m-cycle where the subscripts k n and s n are determined by 0 k n s n and t k n t n < t k n +1 , t s n t n + Φ(α(x n ,t n )) < t s n +1 .
Note that t i+m = t i + T Σ for any i = 0, 1, . . . and where N 1 is a integer. In order to simplify notation, without loss of generality, assume t 0 = 0, N = 1 and denote τ =t 0 − t k 0 .
is completely defined by the two sequencesx 0 , . . . ,x p−1 ,T 0 , . . . ,T p−2 and τ , whose values can be found by solving of the following system of transcendental equations: for n = 1, . . . , p − 1, where Θ j = t j − t k 0 = Remark 1 For a given m-cycle C m (x(t), t n ) of the plant, bounds on the least period p of the corresponding observer cycle C p (x(t),t n ) can be derived from Eq. (2) as below where the number r denotes the least integer greater than r and r is the largest integer not greater than r (i.e., the integer part, or floor, of r ).
Then, one concludes that, for each n = 1, . . . , m−1, it applies that r n = e D(t n −t n−1 ) r n−1 −λ n−1 B , Therefore, for each n = 1, . . . , m, it follows that Considering only cycles of the same periodicity in the observer as in the plant results in the following simplification.
Proof Note that for n = m one has n−1 i=0T i . As a result, the sequencex 0 , . . . ,x m−1 can be reconstructed as follows: The situation covered by Corollary 1 is, for instance, encountered further in the paper and illustrated in Fig. 3a.
Given the results of [10], the scenarios of 1-cycle and 2-cycle featuring one or two GnRH pulses fired by the plant feedback in the least period of the solution are particularly relevant to the case of Te regulation and considered below.
1-cycle, p = m = 1: The Eq. (13) take the following form: where The matrix M is always non-singular for a Hurwitz D but, as soon as det D = 0 and the continuous observer dynamics are marginally stable, the inverse does not exist. Notably, non-singularity of M can always be secured by a proper choice of K since (A, C) is observable. 2-cycle, p = m = 2: Equation (13) are reduced to: where T Σ = T 0 + T 1 .

Stability of periodic observer solutions
Stability conditions for an m-cycle C m (x(t),t n ) of observer Eqs. (4)- (6) satisfying can be obtained through analyzing the Jacobian of the pointwise map iterated m times F m = Q (m) evaluated at the fixed point of map Eq. (11). Then, the m-cycle is locally asymptotically stable if the Jacobian F m = Q (m) is Schur stable [15,17]. The calculation of the Jacobian can be simplified by applying the chain rule, as

Theorem 2
The Jacobian of Q at the point (x n ,t n ) is calculated as Here,Φ n ,F n ,Φ n andF n denote Φ(α(x n ,t n )), F(x n ) and their derivatives, respectively, and x + i = x(t + i ). The indices k and s obey the inequality t k t n < t k+1 , t s t n +Φ n < t s+1 .
Proof By direct calculations along the lines of the proof of Theorem 2 in [15].
If (x(t),t n ) is a synchronous mode, i.e., (x(t),t n ) = (x(t), t n ), then the Jacobian of Q at the point (x n , t n ) is calculated as where (J n ) 11 = Φ (Cx n )Ax n+1 R + e DΦ(Cx n ) I + F (Cx n )BC , Local asymptotic stability implies that if the point (x n ,t n ) is in a sufficiently small neighborhood of the point (x n , t n ) then, after a firings of the observer discrete feedback, the point (x n+a ,t n+a ) will be in a smaller neighborhood of the point (x n+a , t n+a ), for some a 1. Yet, in [22], it was shown that the basin of attraction of a synchronous mode in the vicinity of the point (x n , t n ) is considerably asymmetrical relative to the jump instant t n . Even small overestimates of the plant feedback firing time (i.e.,t n > t n ) interrupt and postpone the convergence. The term "overjump" will be used to denote a situation when the estimated firing timet n+1 lies in a right neighborhood of t n+1 , whilet n belongs to a left neighborhood of t n . Since the function Φ(·) is increasing, the overjumps are generally more likely to occur for high values of k d . This may lead to the appearance of undesired high-periodic or chaotic modes of the observer that can be revealed and investigated by bifurcation analysis; see Sect. 7.
After finding a feasible set of the gains K and k d for which the synchronous mode is attractive in the whole D, one has to determine the suitable values of K and k d ensuring the highest achievable convergence rate. Since the greatest challenge for the observation problem at hand is the unknown initial firing time t 0 (and therefore all the following instants t n ) of the system, the observer convergence rate can be characterized by the first time instant whent n comes into ε f -neighborhood of t n and never leaves it, i.e.,

Bifurcation analysis
The normal stationary operation regime for the observer is that of synchronous mode. All other modes have to be excluded by design, while the convergence to a synchronous mode has to be guaranteed from any point of the plant initial condition set D at, preferably, the fastest achievable rate.
In the present section, the nonlinear dynamics of the observer at hand are investigated specifically with respect to variation in the observer gain k d . As the analysis below demonstrates, there is a well-defined interval of the gain values that correspond to a synchronous mode of the observer with sufficient basin of attraction and convergence rate. The existence of this interval is guaranteed by the fact that the hybrid state estimation error asymptotically vanishes for a sufficiently large value of the continuous gain K with k d = 0 if the following condition holds (see Theorem 4 in [15]): The necessary condition above reminds once again that the design of observer Eqs. (4)-(6) is performed with respect of to a certain solution of the plant rather than with respect to the model. It also hints to the conclusion that plants with a steep frequency modulation function are more difficult to deal with. Notice also that condition Eq. (18)  In order to discern the nonlinear dynamical phenomena in the observer by means of bifurcation analysis, assume the following parameter values: For the selected parameter values, plant Eq. (3) has a stable 2-cycle depicted in Fig. 2 and inequality Eq. (18) holds. Figure 3a and b demonstrates an example of the stable asynchronous mode arising for k d = −72. The figures depict firing mismatch in the discrete state variable (Fig. 3a) as well as nonzero error in the estimates of the continuous plant states (Fig. 3b).
The hybrid observer dynamics are described by Eqs. (4)- (6). Consider the gain matrix K of the following form: where k c > 0 and renders the matrix D Hurwitz stable. The motivation for restricting the structure of K . Note that since the plant and the observer each have two impulses in the least period, the observation error has four impulses in its least period and the performance consequences of it are discussed in detail in [17]. Here, it suffices to note that, under this restriction, the observer demonstrates good performance [15,17]. The guidelines for selecting k c are given in Sect. 8. Here, k c = 1 is assumed. The choice of the discrete gain k d directly impacts the discrete dynamics of the observer and is, therefore, more significant.

Discrete gain
In order to discern the mechanisms of the transition from asynchronous to synchronous mode, select k d as the bifurcation parameter. Figure 4 displays the results of one-dimensional bifurcation analysis for −150.0 < k d < 100.0. Figure 4a shows the bifurcation diagram obtained through direct simulation and by means of continuation techniques using transcendental Eq. (15). These equations allow to determine not only the stable cycles but as well the unstable ones. The two horizontal lines denoted by number 1 in Fig. 4a represent the observed 2-cycle of the plant. Figure 4b, c illustrates the variation of multipliers of the asynchronous and synchronous 2-cycles.
For the values of the gain factor k d < k sn d , mapping Eq. (11) shows a pair of asynchronous 2-cycles (Fig. 4a), one of which is stable (drawn in solid lines and denoted by number 3), while the other, drawn in dashed lines and denoted by number 2, is a saddle one.
As the gain k d increases, the stable 2-cycle (denoted by number 3) merges with the saddle 2-cycle (2) and disappears in a saddle-node bifurcation at the point k sn d (Fig. 4a). The variation of the multipliers ρ 1,2 and ρ u 1 for the stable and saddle asynchronous 2-cycles, respectively, is shown in Fig. 4b, c. To better illustrate the later transition, a magnified section of Fig. 4a outlined by the dashed rectangle is presented in Fig. 4c. As illustrated in this figure, with the increasing gain factor k d , the stable focus 2-cycle transforms into the stable node 2-cycle, as a pair of complex-conjugated multipliers ρ 1,2 = μ ± iω for the stable 2-cycle become real ρ 1 and ρ 2 . At the saddle-node bifurcation point k sn d , the largest (in the absolute value) multipliers ρ 1 and ρ u 1 of the stable and saddle 2-cycles become equal to + 1.
When the gain factor k d passes the value k d = k sn d , an abrupt transition to chaos takes place, as evidenced by Fig. 4a, c, and the observer does not exhibit stable periodic dynamics in the interval between the points k sn d and k sn * d . In the region k sn d < k d < k sn * d , the largest Lyapunov exponent Λ 1 becomes positive, signaling the development of chaotic dynamics (Fig. 4d). The second Lyapunov exponent Λ 2 , on the other hand, is negative everywhere in the considered interval −150.0 < k d < 100.0.
With the further increase in k d , the system enters the 2-cycle window. The saddle-node bifurcation at the left edge k sn * d of the region k sn d < k d < k sn * d produces the saddle and stable asynchronous 2-cycles, drawn in dashed line and denoted by number 3, and drawn in solid line and denoted by number 2, respectively. When crossing the point k * d with increasing k d , the stable asynchronous 2-cycle merges with the saddle synchronous 2-cycle in a "transcritical"-like bifurcation: at the point k * d , one of the multipliers ρ * 1 and ρ * * 1 of the saddle synchronous and stable asynchronous 2-cycles, respectively, becomes equal to + 1 (Fig. 4b, c). After this bifurcation, the saddle synchronous 2-cycle transforms into a stable one (denoted by number 4). The variation of the multipliers for the synchronous (ρ * 1,2 ) and asynchronous (ρ * * 1,2 ) 2-cycles is shown in Fig. 4b. The considered bifurcation is similar to the classical transcritical bifurcation. However, dissimilar to the transcritical bifurcation, this transition does not involve the appearance of an unstable asynchronous 2-cycle.
As illustrated in Fig. 4b, when crossing the point k c * d with increasing gain parameter k d , the stable node synchronous 2-cycle turns into stable focus (the real multipliers ρ * 1,2 become complex ρ * 1,2 = μ * ± iω * ). The domain that falls to the right of the point k h d is a region of bistability, where the stable synchronous 2-cycle coexists with chaotic or high-periodic attractors. In the domain of chaotic dynamics, there exists a usual dense set of periodic windows (Fig. 4a). Figure 4e depicts the magnified part of the diagram presented in Fig. 4d of the largest Lyapunov exponent Λ 1 in the region k d > k h d . As illustrated in Fig. 4e, the largest Lyapunov exponent Λ 1 is positive in most of the considered interval, and this Lyapunov exponent becomes negative in the windows of periodicity.
It can be concluded that the basins of attraction of the coexisting motions are delineated by the stable manifold of the saddle asynchronous 2-cycle denoted by number 2 in Fig. 4a, which appears at the point k sn d in a saddle-node bifurcation. Figure 4f displays the magnified part of the bifurcation diagram that is outlined by the rectangle in Fig. 4a. Inspection of Fig. 4a, f reveals that, as the parameter k d increases, the saddle asynchronous 2-cycle approaches the stable synchronous 2-cycle, when of the gain factor k d changes the sign.
In the region of bistability k d > k h d , the stable synchronous and saddle asynchronous 2-cycles are very close to each other. Here, the system may choose different motions depending on initial conditions. In particular, if a certain ratio between the radii of basins of attractions and the magnitude of random disturbances takes place, then either chaotic or high-periodic dynamics occurs. It also becomes possible for the system to exhibit a hard transition from a stable synchronous 2cycle to a chaotic or high-periodic attractor and vice versa. Fig. 6 Variation of the absolute value of the largest multipliers for the synchronous 2-cycle. As illustrated in Fig. 4b, when passing through the value k c * d , the largest (in the absolute value) real multiplier of the 2-cycle becomes equal to another one, which gives rise to a pair of complex-conjugated multipliers Figure 5a depicts a three-dimensional phase space projection of the coexisting stable synchronous 2-cycle (denoted by C 2 ) and stable 31-cycle (denoted by C 31 ) for k d = 84.5. Figure 5b illustrates a two-dimensional projection of the basins of attraction for coexisting 2and 31-cycles. Here, B 2 and B 31 denote the basins of attraction for 2-and 31-cycles, respectively.

Basin of attraction for the synchronous mode
As mentioned earlier, it is important to guarantee the observer convergence to the synchronous mode from all admissible initial estimates, i.e., in the whole set D. For the numerical example at hand, according to Eq. As illustrated in Figs. 4a, b and 5, the synchronous mode is locally asymptotically stable for k d k * d . Figure 6 shows the absolute value of the largest multipliers ρ for the synchronous 2-cycle. This value will be referred to as the spectral radius ρ(F 2 ) of the Jacobian F 2 .
In Fig. 6, one can see that the synchronous mode is still locally stable for k d > k h d since ρ(F 2 ) < 1. However, these values of k d are not suitable for the observer since the basin of attraction for the synchronous mode is too small, and the synchronous 2-cycle either coexists with stable high-periodic asynchronous modes or with a chaotic attractor (Figs. 4 and 5).

Convergence rate to the synchronous mode
When it is guaranteed that the observer estimates converge to the synchronous mode from all the points of D for k * d k d k h d , it is reasonable to select the gain value that ensures fastest observer convergence. The greatest challenge for the hybrid observation problem at hand is the unknown initial firing time t 0 . Then, a suitable gain k d from the interval k * d , k h d should be chosen to minimize the settling time P(ε f ) defined in Eq. (17), for all initial timeŝ t 0 from T = {θ ∈ R : 0 θ T Σ }. From Fig. 7, it can be concluded that for the fastest convergence of the observer to the synchronous mode from all possible initial instantst 0 one should select a value of k d close to the upper bound of the feasible interval, i.e., k h d . Since initial values of the unmeasurable components of the system state x cannot be related tot 0 , one has to select them independently oft 0 , e.g.,x 0 = x 1,0 ,x 2,0 ,x 3,0 . Notice that the initial condition for the measurable component estimatex 3,0 =ẑ(t 0 ) has the greatest impact on the convergence, since it is the modulation parameter and defines the value ofT 0 .

Observer design algorithm
In this section, an algorithm for finding the values of the observer gains K and k d rendering the synchronous mode attractive for all the points of D with an appropriate convergence rate is presented.
Let the structure of K be defined by Eq. (19). First of all, notice that local stability of the synchronous mode is a necessary condition for attractivity of D. Since Eq. (18) holds, for k d = 0, one can find a value k max such that, for any k c > k max , the spectral radius ρ of the Jacobian F 2 calculated at the points of the synchronous mode differs no more than ε > 0 from ρ(F 2 ) for k c = k max . Figure 9 shows the dependence of ρ(F 2 ) on the gains k c and k d . For a fixed k c (0 < k c < k max ), denote by K(k c ) the interval of k d for which the synchronous mode is locally asymptotically stable (ρ(F 2 ) < 1). For example, in Fig. 9, the intervals K(0.5) = (−14.9, 90.1) and K(1) = (−49. 5, 288.7) are highlighted by the red horizontal line segments. One can see that K(k c ) grows with increasing k c . Moreover, since inequality Eq. (18) holds, the spectral radius ρ(F 2 ) → ρ const < 1 as k c → ∞, for any fixed k d . However, the observer convergence rate is low for large values of k c , whereas one can see from  Ift 0 is far from t 0 , then the observer may not approach the synchronous mode at all, e.g., for small k c , the interval of suitable k d becomes very small and settling time increases (Fig. 10).
Denote by K 2 (k c ) ⊆ K(k c ) the interval k d for which the whole subspace of feasible initial conditions D is within the basin of attraction of the synchronous mode. The interval K 2 (k c ) can be found by bifurcation analysis. For k c = 1, the interval is evaluated to K 2 (k c ) = (−49.5, 38.5) (Fig. 4). Notice that it can be sufficient to determine the interval K 3 (k c ) ⊆ K 2 (k c ) that consists of those k d that satisfy P(ε f ) < P * for all initial conditions in D, where P * is some large number. Figure 10 shows the dependence of P(ε f ) on k c an k d for ε f = 1 andt 0 = 120, where the two red horizontal line segments are the intervals of k d such that P(ε f ) < 8000 for k c = 1 and k c = 2, respectively. For k c = 1, it was found that K 3 (1) = (−25.7, 38.3). From Fig. 10, one can conclude that, for large k c , the interval of k d , for which the observer converges to the synchronous mode, increases. However, the interval K 3 (k c ) is bounded, since overjumps may occur for some values oft 0 and large k d . Figure 10 shows that, for each large k c , one can choose k d ensuring fast convergence (i.e., the values in the dark blue area). Yet, with a large k d , convergence cannot be achieved for all possible initialt 0 ∈ T due to overjumps. Therefore, it is reasonable to choose some moderate value of k c that yields a sufficiently large interval of suitable k d such that the synchronous mode is attractive for allt 0 ∈ T . Finally, after finding such a value of k c , a suitable k d that minimizes P(ε f ) overt 0 ∈ T can be chosen from Recapitulating, the following bifurcation analysisbased iterative algorithm for finding the gains K and k d for a given periodic mode in the plant is proposed. For non-periodic (chaotic, quasi-periodic) modes in the plant, the design procedure will be the same, but one has to use the Lyapunov exponents instead of the spectral radius of the Jacobian.
For k d = 0, find the value k max . do: Starting with k c = k max and decreasing k c with some step Δ do the following (for each fixed k c ) Step 1: determine the interval of k d for which the synchronous mode is locally stable, i.e., the spectral radius of the Jacobian is less than one (or the largest Lyapunov exponent is negative) and denote it K; Step 2: exclude the values of k d for which the stable synchronous mode, according to bifurcation analysis, coexists with other stable modes, as a result obtain the set K 1 ⊆ K; Step 3: find a subset K 2 ⊆ K 1 such that the whole subspace of feasible initial conditions D is within the basin of attraction of the synchronous mode; Step 4: find the interval K 3 ⊆ K 2 containing the values of k d for which the settling time P(ε f ) < P * for all initial conditions from D; Step 5: finally, select k d from K 3 that minimizes the criterion P(ε f ).
As a suitable gaink c choose the one, for which the value of P(ε f ) found in Step 5 and averaged over all t 0 ∈ T is minimal. The corresponding value of the gain k d should be taken as suitable k d .

Design example
The main intended application area for the developed observer is the processing of endocrine data where it is supposed to complement or substitute the deconvolution techniques. Naturally, to be implemented, the observer demands a complete plant model. Such individualized models have been obtained from experimental data in [10]. At that stage, no consistent theory for design of the observer was available and pure simulation of the estimated model was performed to demonstrate the fidelity. System identification is not considered in this work, and therefore, the observer is not run on biological data but only on synthetic such.
For the example considered in Sect. 7, one obtains, with ε = 0.02, that k max = 3. Decreasing k c with the step Δ = 0.1 (i.e., N = 30), the calculations specified by Step 1-Step 5 of the proposed observer design algorithm are performed, for each value of k c . As a result, one obtains k c = 1, k d = 38.3 as suitable observer gains.
To illustrate Step 1-Step 5 of the algorithm in Sect. 7, the results of one cycle of it for k c = 1 are provided below: Step 1: the interval K = (−49.5, 288.7); see Figs. 6 and 9; Step 2: the interval K 1 = (−49.5, k h d ) = (−49.5, 38.5); see Fig. 4a; Step 3: an analysis of the basin of attraction of the synchronous mode in the four-dimensional initial conditions space reveals that it includes all the points within D for the interval K 2 = K 1 ; Step 4: for P * = 8000, one can conclude from Fig. 10 that the interval K 3 = (−25.7, 38.3); Step 5: the averaged overt 0 ∈ T minimum value of P(ε f ) for ε f = 1 is achieved for k d = 38.2; see Fig. 7.
The value of P(ε f ) for k c = 1 obtained on Step 5 is minimal for all 0 < k c 3, so k c = 1 and k d = 38.2 are selected as suitable gain values. Figure 11 contains time-domain illustrations of the designed observer performance starting from different initial conditions. Note that the convergence rate of the designed observer is not much dependent on the initial conditions. In fact, this is also seen in Fig. 10, where the rate of convergence practically does not change depending on the components of the statex 0 , since large initial deviations are quickly compensated by an exponential decrease in the first intervals of continuity. Therefore, the initial state ofx does not really matter. The rate of convergence is affected mostly by the discrete componentt 0 , namely the distance to t 0 (especially, for cycles of low periodicity, e.g., 1-cycles and 2-cycles, here the influence of k d is preeminent) and the occurrence of a jump after the first pulse.

Conclusions
A systematic design procedure for a hybrid observer reconstructing the state variables of the impulsive Goodwin's oscillator is proposed. It is based on the numerical bifurcation analysis of a discrete mapping capturing the observer state transition from one firing of the cumulative impulsive sequence of the observer and plant feedback to another one. The design procedure produces the values of the observer continuous and discrete gains that guarantee the asymptotic convergence of the hybrid state estimation error to zero from all feasible initial conditions at a highest possible rate.
The problem of evaluating periodic solutions in the observer is reduced to solving a system of analytically derived transcendental equations with respect to the points of the periodic orbit of the discrete mapping. This allows to not only evaluate the stable cycles but as well the unstable ones.
A detailed investigation of the bifurcation phenomena in transition from asynchronous to synchronous observer mode with respect to the discrete gain is presented. In particular, the transition from a stable asynchronous 2−cycle to a stable synchronous 2-cycle can occur through a "transcritical"-like bifurcation in which the saddle synchronous 2-cycle transforms into a stable one. However, unlike to the transcritical bifurcation, this transition does not involve the appearance of an unstable asynchronous 2-cycle. For higher values of gain, the system displays a region of bistability, where the stable synchronous 2-cycle coexists with chaotic or high-periodic attractors. Since the stable coexisting synchronous and saddle asynchronous 2-cycles are very close to each other, the observer is sensitive to noise this region. An arbitrarily small level of noise may cause a sudden and unpredictable transition between the attractors.
The design procedure is illustrated by a simulation example that confirms efficacy of the observer.