Switching Polynomial Systems

• SPSs are hybrid systems that switch among finitely many polynomial modes, integrating continuous dynamics with discrete decision rules.
• They enable rigorous stability analysis using common Lyapunov functions and SOS certificates, while supporting optimal switching design and safe controller synthesis.
• Applications include invariant set computation, bifurcation analysis, and system identification, leveraging computational methods from LP relaxations to formal proof techniques.

Switching Polynomial Systems (SPSs) are hybrid dynamical systems whose evolution switches among finitely many polynomial modes. In the standard continuous-time form, one writes

$\dot{x}(t)=f_{\sigma(t)}(x(t)),$

with polynomial vector fields $f_i:\mathbb{R}^n\to\mathbb{R}^n$ and a switching signal $\sigma$; a discrete-time analogue is $x_{k+1}=f_{\sigma(k)}(x_k)$. In planar piecewise formulations, switching occurs across a manifold such as $x=0$, $y=0$, or $\mathbb{S}^1$, whereas in state-dependent formulations the active mode is determined by polynomial sign conditions. A more specialized construction replaces the switching signal itself by polynomial weights, notably Bernstein interpolants, so that the hybrid dynamics are recast as a single polynomial vector field with nonnegative convex coefficients. Across stability theory, controller synthesis, optimal switching, invariant-set computation, formal verification, bifurcation analysis, and identification from data, SPSs combine the algebraic tractability of polynomial systems with the combinatorial structure of switching (Michel, 2014, Ahmadi et al., 2017).

1. Core definitions and modeling variants

A common definition of an SPS is a finite family of polynomial vector fields $\{f_i(x)\}_{i=1}^m$ together with a switching signal $\sigma(t)\in\{1,\ldots,m\}$, typically piecewise constant, so that

$\dot{x}(t)=f_{\sigma(t)}(x(t)).$

In the linear special case, $f_i:\mathbb{R}^n\to\mathbb{R}^n$0, yielding switched linear systems; homogeneous polynomial vector fields and planar polynomial vector fields are treated as important subclasses in the Lyapunov and SOS literature (Ahmadi et al., 2017).

A second modeling line treats SPSs as semi-algebraic hybrid systems. Here the discrete modes $f_i:\mathbb{R}^n\to\mathbb{R}^n$1 are equipped with polynomial flows $f_i:\mathbb{R}^n\to\mathbb{R}^n$2, mode domains

$f_i:\mathbb{R}^n\to\mathbb{R}^n$3

and guards

$f_i:\mathbb{R}^n\to\mathbb{R}^n$4

Under the standing assumptions of polynomial data, complete flows, and identity resets, switching-controller synthesis can be formulated as the refinement of domains and guards within the original hybrid automaton (Kapur et al., 2013).

A third modeling line emphasizes state-dependent switching surfaces. In continuous-time identification settings, SPSs are written as

$f_i:\mathbb{R}^n\to\mathbb{R}^n$5

with polynomial modes $f_i:\mathbb{R}^n\to\mathbb{R}^n$6, where $f_i:\mathbb{R}^n\to\mathbb{R}^n$7 is a polynomial feature map of monomials up to degree $f_i:\mathbb{R}^n\to\mathbb{R}^n$8. The switching rule is encoded by polynomial surfaces

$f_i:\mathbb{R}^n\to\mathbb{R}^n$9

whose signs partition the state space into mode regions $\sigma$0 (Iwasaki et al., 29 Sep 2025).

In planar nonsmooth dynamics, SPSs are often piecewise-defined polynomial vector fields separated by a switching manifold. For two-zone systems,

$\sigma$1

with Filippov conventions on the switching set. This framework includes switching Liénard systems, $\sigma$2-equivariant cubic switching systems, and piecewise polynomial holomorphic or antiholomorphic systems with straight or circular switching manifolds (Chen et al., 2023, Chen et al., 2024, Rondón et al., 28 Apr 2026).

A more specialized usage arises in the Bernstein-based formulation of “B-switched systems,” where switching itself is polynomially encoded: $\sigma$3 Here the switching signal is replaced by Bernstein-polynomial interpolants, producing a unified polynomial dynamics with convex weights (Michel, 2014).

2. Lyapunov theory, SOS certificates, and polynomial stability criteria

The central stability question for SPSs is the existence of a Lyapunov certificate compatible with all modes. For a common Lyapunov function $\sigma$4, the continuous-time inequalities take the form

$\sigma$5

with $\sigma$6 on the punctured domain. In discrete time, one analogously requires

$\sigma$7

When such a common $\sigma$8 exists, stability holds under arbitrary switching; multiple Lyapunov functions with dwell-time or average dwell-time constraints constitute an alternative when a common certificate is unavailable (Michel, 2014).

For several SPS subclasses, converse SOS results are available. For homogeneous polynomial vector fields, the existence of a homogeneous polynomial Lyapunov function $\sigma$9 with $x_{k+1}=f_{\sigma(k)}(x_k)$0 positive definite implies the existence of another homogeneous polynomial Lyapunov function $x_{k+1}=f_{\sigma(k)}(x_k)$1 such that $x_{k+1}=f_{\sigma(k)}(x_k)$2 is SOS and $x_{k+1}=f_{\sigma(k)}(x_k)$3 is SOS. For planar polynomial vector fields, the same conclusion holds under the mild condition that the top homogeneous component $x_{k+1}=f_{\sigma(k)}(x_k)$4 has no zeros. For stable switched linear systems, the result sharpens to an if-and-only-if statement: $x_{k+1}=f_{\sigma(k)}(x_k)$5 is asymptotically stable under arbitrary switching if and only if there exists a common homogeneous polynomial Lyapunov function $x_{k+1}=f_{\sigma(k)}(x_k)$6 such that $x_{k+1}=f_{\sigma(k)}(x_k)$7 is SOS and, for all modes, $x_{k+1}=f_{\sigma(k)}(x_k)$8 is SOS, all positive definite (Ahmadi et al., 2017).

The same work also isolates a useful weakening of the usual SOS requirement. If $x_{k+1}=f_{\sigma(k)}(x_k)$9 is SOS and $x=0$0 is SOS, then global asymptotic stability follows; moreover, “$x=0$1 SOS” is strictly weaker than “$x=0$2 SOS,” and in the planar case it is lossless in conservatism because nonnegative bivariate forms are SOS. This modifies the usual SDP search space and can reduce decision variables and equality constraints (Ahmadi et al., 2017).

A different route to polynomial Lyapunov functions is the Lyapunov Differential Equation hierarchy for switched linear systems. At hierarchy level $x=0$3, the $x=0$4-th Kronecker power of the state evolves linearly under lifted matrices

$x=0$5

and feasibility of

$x=0$6

yields a homogeneous polynomial Lyapunov function of degree $x=0$7. After reduction to a unique-monomial basis $x=0$8, one obtains

$x=0$9

which is computationally competitive with SOS while remaining equivalent to homogeneous SOS Gram representations (Abate et al., 2019).

3. Bernstein convexification and formal verification

The Bernstein-based approach begins from the univariate basis

$y=0$0

and the Bernstein polynomial

$y=0$1

with uniform convergence to $y=0$2 on $y=0$3 for continuous functions. Multivariate extensions are defined on $y=0$4 by tensor-product Bernstein bases. Because the coefficients are nonnegative and sum to one, Bernstein polynomials preserve positivity and remain within the convex hull of their coefficients (Michel, 2014).

These properties are used to replace a switching signal by polynomial weights. For two modes with a state-dependent sign guard, the construction takes

$y=0$5

with the guard argument scaled into $y=0$6. The resulting B-switched system is a single polynomial vector field with nonnegative convex coefficients. In the numerical illustrations, low-degree Lyapunov functions are combined with switching interpolants of degree $y=0$7, and the method is analyzed on compact hypercubes scaled to $y=0$8 (Michel, 2014).

This polynomial encoding supports formal verification in HOL Light. The workflow is: choose $y=0$9, compute $\mathbb{S}^1$0 for the Bernstein-convexified dynamics, bound $\mathbb{S}^1$1 either through interval/Taylor methods or Bernstein coefficients, and invoke HOL Light’s verify_ineq on the domain box. The paper’s illustrative theorem-proving example certifies $\mathbb{S}^1$2 on $\mathbb{S}^1$3 for a Hurwitz linear system with

$\mathbb{S}^1$4

The same template is proposed for SPSs once switching has been polynomially encoded (Michel, 2014).

The method’s scope is explicit. Modes must be polynomial or polynomially approximated on the chosen domain; discontinuous switching signals such as exact sign functions are only approximated; the analysis is local to compact boxes; and higher dimensions or higher degrees increase both coefficient counts and proof complexity. The trade-off is a clean path from polynomial encoding to machine-checked stability certificates (Michel, 2014).

4. Optimal switching design and switching-controller synthesis

For optimal control, SPSs can be treated as hybrid systems whose switching sequence is the design variable. The occupation-measure formulation introduces, for each mode $\mathbb{S}^1$5, a nonnegative measure $\mathbb{S}^1$6 on $\mathbb{S}^1$7, together with initial and terminal measures $\mathbb{S}^1$8. The dynamics are encoded weakly by the Liouville equation

$\mathbb{S}^1$9

and the time-marginal normalization enforces that the sum of mode occupations has Lebesgue time as marginal. This yields an infinite-dimensional linear program whose moment-SOS hierarchy provides a nondecreasing sequence of lower bounds $\{f_i(x)\}_{i=1}^m$0 under compactness and Archimedean assumptions (Henrion et al., 2013).

The finite-dimensional SDP relaxation is built from moment and localizing matrices, with support constraints for semialgebraic state sets and time intervals. From the resulting moments, one can extract almost optimal switching schedules by reconstructing time marginals, differentiating them to reveal jump atoms, and recovering switching times by Prony/ESPRIT-type procedures or generalized eigenvalue decompositions. The framework is convex, provides certified global lower bounds, and accommodates chattering in the relaxed problem (Henrion et al., 2013).

A separate synthesis tradition constructs switching controllers by continuous invariant generation. Given a hybrid automaton $\{f_i(x)\}_{i=1}^m$1 and safety sets $\{f_i(x)\}_{i=1}^m$2, the goal is to synthesize a refinement $\{f_i(x)\}_{i=1}^m$3 such that $\{f_i(x)\}_{i=1}^m$4, all trajectories remain safe, and the system is non-blocking. The key theorem states that if, for each mode $\{f_i(x)\}_{i=1}^m$5, $\{f_i(x)\}_{i=1}^m$6 is a continuous invariant of $\{f_i(x)\}_{i=1}^m$7, where

$\{f_i(x)\}_{i=1}^m$8

then the refined automaton solves the synthesis problem (Kapur et al., 2013).

The synthesis constraints can be solved symbolically by quantifier elimination, yielding relative completeness with respect to a chosen template family, or numerically by SOS relaxations and template polyhedra. The SOS route encodes invariance and safety inclusions through polynomial multipliers and semidefinite constraints; the template-polyhedra route specializes to linear dynamics and convex invariants, using generalized Farkas conditions that lead to bilinear matrix inequalities. The paper’s nuclear reactor example exhibits all three routes—QE, SOS, and template polyhedra—within a single SPS synthesis problem (Kapur et al., 2013).

5. Invariant sets, constraints, and domains of attraction

A closely related problem is the computation of admissible invariant sets under polynomial state constraints. For discrete-time switching linear systems

$\{f_i(x)\}_{i=1}^m$9

with arbitrary switching and constraint set

$\sigma(t)\in\{1,\ldots,m\}$0

Veronese lifting maps the polynomial inequalities to linear inequalities in lifted coordinates $\sigma(t)\in\{1,\ldots,m\}$1, while the lifted dynamics remain linear: $\sigma(t)\in\{1,\ldots,m\}$2 This converts a non-convex semi-algebraic constraint problem into a polyhedral invariance problem in lifted space (Athanasopoulos et al., 2015).

Under the assumptions that the joint spectral radius satisfies $\sigma(t)\in\{1,\ldots,m\}$3, that $\sigma(t)\in\{1,\ldots,m\}$4 is compact, and that the origin lies in its interior, the maximal admissible positively invariant set coincides with the constrained domain of attraction. The robust predecessor operator

$\sigma(t)\in\{1,\ldots,m\}$5

supports backward-reachability iterations that converge in finitely many steps. Three exact variants are given, depending on whether one initializes from $\sigma(t)\in\{1,\ldots,m\}$6, from $\sigma(t)\in\{1,\ldots,m\}$7 with a bounding set, or from a bounded polyhedral initialization (Athanasopoulos et al., 2015).

The algorithmic core is LP-based: predecessors of lifted polyhedra remain polyhedral, redundant inequalities can be removed by linear programs, and optional SOS-based SDP post-processing can simplify the lowered semi-algebraic description. The same paper also discusses implications for genuine SPSs with polynomial mode maps: exact finite-dimensional closure generally fails, so one is led either to Carleman-type liftings with truncation or to direct SOS certification of polynomial predecessor relations. In that sense, exact Veronese lifting is a special feature of linear modes, whereas its extension to SPSs becomes approximate or relaxation-based (Athanasopoulos et al., 2015).

6. Nilpotent centers, bifurcations, and limit cycles

A substantial branch of SPS theory concerns planar oscillatory dynamics. For cubic switching polynomial Liénard systems split by $\sigma(t)\in\{1,\ldots,m\}$8,

$\sigma(t)\in\{1,\ldots,m\}$9

a higher-order Poincaré–Lyapunov method is developed to handle nilpotent singularities, where classical linear-center analysis is unavailable. The paper derives eight explicit nilpotent-center conditions for cubic switching Liénard systems, characterizes global nilpotent centers in two cases, and proves that under one center condition proper Bogdanov–Takens-type perturbations yield at least five small-amplitude crossing limit cycles, a new lower bound of cyclicity for that class (Chen et al., 2023).

For planar $\dot{x}(t)=f_{\sigma(t)}(x(t)).$0-equivariant cubic switching systems separated by the $\dot{x}(t)=f_{\sigma(t)}(x(t)).$1-axis, with two symmetric nilpotent singular points at $\dot{x}(t)=f_{\sigma(t)}(x(t)).$2, a generalized e-Lyapunov method is used to derive six explicit bi-center conditions. Building on those conditions and combining generalized Lyapunov constants with Bogdanov–Takens bifurcation theory and a symmetric pseudo-Hopf perturbation, the paper proves the existence of at least $\dot{x}(t)=f_{\sigma(t)}(x(t)).$3 small-amplitude limit cycles around the nilpotent bi-center, increasing the lower bound of cyclicity from $\dot{x}(t)=f_{\sigma(t)}(x(t)).$4 to $\dot{x}(t)=f_{\sigma(t)}(x(t)).$5 for that cubic switching class (Chen et al., 2024).

A geometrically different setting arises when the switching manifold is the unit circle $\dot{x}(t)=f_{\sigma(t)}(x(t)).$6. Möbius transformations are used to map circular discontinuities to straight-line discontinuities while preserving periods, stability, and algebraic structure. For piecewise polynomial holomorphic systems, second-order averaging yields lower bounds

$\dot{x}(t)=f_{\sigma(t)}(x(t)).$7

and low-degree weak-focus constructions yield $\dot{x}(t)=f_{\sigma(t)}(x(t)).$8 for $\dot{x}(t)=f_{\sigma(t)}(x(t)).$9 with $f_i:\mathbb{R}^n\to\mathbb{R}^n$00. For piecewise polynomial antiholomorphic systems, Hamiltonian and resultant arguments give global upper bounds: at most $f_i:\mathbb{R}^n\to\mathbb{R}^n$01 limit cycles in the linear case and $f_i:\mathbb{R}^n\to\mathbb{R}^n$02 in the quadratic case. The same work proves a rigidity theorem excluding crossing limit cycles when both components admit classical holomorphic normal forms at the origin, and it constructs explicit algebraic limit cycles in the circular setting (Rondón et al., 28 Apr 2026).

Together, these results show that SPSs are not only a framework for stability certification but also a setting in which switching geometry, symmetry, nilpotent degeneracy, and Filippov phenomena materially alter center conditions and cyclicity bounds.

7. Identification from data and computational developments

Recent work addresses SPSs as objects of system identification. In the state-dependent formulation,

$f_i:\mathbb{R}^n\to\mathbb{R}^n$03

the joint estimation of mode assignments and polynomial mode dynamics is first posed as a mixed-integer program with one-hot variables $f_i:\mathbb{R}^n\to\mathbb{R}^n$04. To address scalability, the paper develops both semidefinite moment relaxations and simplex relaxations. The moment approach admits a flat-extension tightness criterion—rank-one moment matrices recover exact one-hot assignments—while the simplex relaxation is shown to be generically one-hot under mild nondegeneracy and absolute continuity assumptions (Iwasaki et al., 29 Sep 2025).

The computational scheme is bilevel and alternating. With $f_i:\mathbb{R}^n\to\mathbb{R}^n$05 fixed, mode assignment is solved by an order-1 moment SDP or, more cheaply, by a simplex LP under an $f_i:\mathbb{R}^n\to\mathbb{R}^n$06 residual model. With assignments fixed, the polynomial coefficients $f_i:\mathbb{R}^n\to\mathbb{R}^n$07 are updated by LP under box constraints. The resulting iterate sequence has monotone decrease of the $f_i:\mathbb{R}^n\to\mathbb{R}^n$08 objective, and any accumulation point is blockwise optimal. After mode labels are obtained, switching boundaries are reconstructed by a soft-margin polynomial classification LP,

$f_i:\mathbb{R}^n\to\mathbb{R}^n$09

with $f_i:\mathbb{R}^n\to\mathbb{R}^n$10-regularization on the classifier coefficients (Iwasaki et al., 29 Sep 2025).

The reported nonlinear SPS experiment is a switching quartic oscillator with energy threshold $f_i:\mathbb{R}^n\to\mathbb{R}^n$11, polynomial modes of degree up to $f_i:\mathbb{R}^n\to\mathbb{R}^n$12, and batch size $f_i:\mathbb{R}^n\to\mathbb{R}^n$13. Over $f_i:\mathbb{R}^n\to\mathbb{R}^n$14 iterations, the SDP-based mode-assignment stage required $f_i:\mathbb{R}^n\to\mathbb{R}^n$15 s total with median $f_i:\mathbb{R}^n\to\mathbb{R}^n$16 s/iter, whereas the simplex LP required $f_i:\mathbb{R}^n\to\mathbb{R}^n$17 s total with median $f_i:\mathbb{R}^n\to\mathbb{R}^n$18 s/iter. With boundary recovery degree $f_i:\mathbb{R}^n\to\mathbb{R}^n$19, $f_i:\mathbb{R}^n\to\mathbb{R}^n$20, $f_i:\mathbb{R}^n\to\mathbb{R}^n$21, and $f_i:\mathbb{R}^n\to\mathbb{R}^n$22, the reported pointwise performance was velocity RMSE $f_i:\mathbb{R}^n\to\mathbb{R}^n$23, mode accuracy $f_i:\mathbb{R}^n\to\mathbb{R}^n$24, and mIoU $f_i:\mathbb{R}^n\to\mathbb{R}^n$25; on $f_i:\mathbb{R}^n\to\mathbb{R}^n$26 rollouts, RMSE was $f_i:\mathbb{R}^n\to\mathbb{R}^n$27, final error $f_i:\mathbb{R}^n\to\mathbb{R}^n$28, and max error $f_i:\mathbb{R}^n\to\mathbb{R}^n$29 (Iwasaki et al., 29 Sep 2025).

This identification framework complements earlier SPS literature centered on analysis, synthesis, and certification. It retains polynomial structure at the level of both dynamics and switching boundaries, while replacing symbolic or theorem-proving workflows by convex relaxations and data-driven estimation (Iwasaki et al., 29 Sep 2025).