{'fraction_non_alphanumeric': 0.03733145502113101, 'fraction_numerical': 0.00823103240088549, 'mean_word_length': 5.097055214723927, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Model-Based Reinforcement Learning for Stochastic Hybrid Systems

Optimal control of general nonlinear systems is a central challenge in automation. Enabled by powerful function approximators, data-driven approaches to control have recently successfully tackled challenging robotic applications. However, such methods often obscure the structure of dynamics and control behind black-box over-parameterized representations, thus limiting our ability to understand closed-loop behavior. This paper adopts a hybrid-system view of nonlinear modeling and control that lends an explicit hierarchical structure to the problem and breaks down complex dynamics into simpler localized units. We consider a sequence modeling paradigm that captures the temporal structure of the data and derive an expectation-maximization (EM) algorithm that automatically decomposes nonlinear dynamics into stochastic piecewise affine dynamical systems with nonlinear boundaries. Furthermore, we show that these time-series models naturally admit a closed-loop extension that we use to extract local polynomial feedback controllers from nonlinear experts via behavioral cloning. Finally, we introduce a novel hybrid relative entropy policy search (Hb-REPS) technique that incorporates the hierarchical nature of hybrid systems and optimizes a set of time-invariant local feedback controllers derived from a local polynomial approximation of a global state-value function.


I. INTRODUCTION
The class of nonlinear dynamical systems governs a vast range of real-world applications and underpins the most challenging problems in classical control, and reinforcement learning (RL) [1], [2]. Recent developments in learning-for-control have pushed towards deploying more complex and highly sophisticated representations, e.g., (deep) neural networks and Gaussian processes, to capture the structure of both dynamics and controllers. This trend led to unprecedented success in the domain of RL [3], and can be observed in both approximate optimal control [4]- [6] and approximate value and policy iteration algorithms [7]- [9].
However, before the latest successful revival of neural networks in control and robotics applications, research focused on different paradigms for solving difficult control tasks. One interesting concept relied on decomposing nonlinear structures  of dynamics and control into simpler local (linear) components, each responsible for an area of the state-action space. This decomposition is done to preserve interpretability and favorable mathematical properties studied over decades in classical control theory, such as local linear-quadratic assumptions [10]. Instances of this abstraction can be found in the control literature under the labels of hybrid systems or switched models [11]- [14], while in the machine and reinforcement learning communities, the terminology of switching dynamical systems (SDS) and switching state-space models (SSM) is more widely adopted [15]- [18]. Building on this vision, we present in this work a view of data-driven automatic system identification and learning of composite control from the perspective of hybrid systems and switching linear dynamics. We are motivated by recent in-depth analysis of piecewise linear (PWL) activation functions such as rectified linear units (ReLU) [19]- [23], which shows that such representations effectively divide the input space into linear sub-regions. This insight highlights the hierarchical structure hidden between a neural network's input and output layers and supports viewing them as approximators that rely on local experts. We take this interpretation as an impulse to follow up on ideas from optimal control [24], [25], and reinforcement learning [26], [27] that deviate from fully differentiable paradigms and investigate whether simpler, hybrid discrete-continuous representations may be sufficient for solving certain tasks. Furthermore, the interest in hybrid systems as graphical models is motivated by favorable properties inherent in such representations. On the one hand, hybrid systems allow the modeling of discrete events, hard nonlinearities, and regiondependent noise. On the other hand, sequence models carry over the advantages of system identification via Bayesian inference and naturally include built-in time recurrent dynamics, which capture correlations over extended time horizons.
This paper consolidates prior work on a hierarchical decomposition of nonlinear dynamics [28] and introduces a novel reinforcement learning algorithm for optimizing local polynomial policies and value functions. In the upcoming sections, we review the literature and highlight the intersection points between prominent paradigms in control and machine learning w.r.t. hybrid systems. Then, we introduce the notation of stochastic switching models and the infinite horizon hybrid control problem. Next, we derive a maximum a posteriori expectation-maximization (EM) algorithm for inferring the probabilistic hybrid dynamics. We use this inference procedure in three different scenarios. First, to perform automatic decomposition of nonlinear open-loop dynamics into switching linear regimes with arbitrary boundaries. Second, to deconstruct state-of-the-art nonlinear expert controllers into simpler local polynomial policies. Finally, we embed the EM procedure into a hybrid policy search algorithm with an explicit discretecontinuous structure. We use this approach to learn hierarchical piecewise polynomial approximations of global value functions and feedback controllers. We empirically evaluate the learned models and policies on a set of numerical examples of stochastic hybrid and nonlinear systems.

II. RELATED WORK
This section reviews work related to the modeling and control of hybrid systems and highlights connections and parallels between approaches stemming from the control and machine and reinforcement learning literature.
Hybrid systems have been extensively studied in the control community and are widely used in real-world applications [29], [30]. For research on the topic of hybrid system identification, we refer to survey work in [31], and [32]. There, the authors focus on piecewise affine (PWA) systems and introduce taxonomies of different representations and procedures commonly used for identifying sub-regimes of dynamics, ranging from algebraic approaches [33] to mixed-integer optimization [34], and Bayesian methods [35]. Furthermore, hybrid system identification techniques for piecewise nonlinear systems have been developed based on sparse optimization [36] and kernel methods [37]. Finally, it is worth noting that the majority of literature considers deterministic mode-switching events with exceptions in [38], [39].
Research in the area of optimal control for hybrid systems stretches back to the seminal work of [40], which highlights the possibility of general nonlinear control by considering piecewise linear systems. In [41], an overview of control approaches for piecewise affine switching dynamics is presented.
The authors categorize the literature by distinguishing between driven and un-driven systems with externally or internally forced switching mechanisms. Given the global nonlinear behavior of switched systems, the bulk of optimal control approaches in this area focus on nonlinear model predictive control (MPC) [42]. Here we highlight the influential work in [43], which formulates the optimal control problem as a mixed-integer quadratic program (MIQP). This approach was later extended in [44] and [45] to solve multi-parametric MIQP and arrive at time-variant local linear state-feedback controllers and local quadratic value functions with affine boundaries. Recently, more efficient formulations of trajectorycentric hybrid control have been proposed [46], which leverage modern techniques from mixed-integer and disjunctive programming and tackle large-scale problems.
Hybrid representations also play a central role in datadriven, general-purpose process modeling and state estimation [47], [48], where different classes of stochastic hybrid systems serve as powerful generative models for complex dynamical behaviors [49]- [51]. The dominant paradigm in this domain has been that of probabilistic graphical models (PGM), more specifically, hybrid dynamic Bayesian networks (HDBN) for temporal modeling [52], [53]. However, one crucial contribution of recent Bayesian interpretations of switching systems is rooted in the Bayesian nonparametric (BNP) view [54]- [57]. This perspective theoretically allows for an infinite number of components, thus dramatically increasing the expressiveness of such models. Given the limited scope of this review section, we highlight only recent contributions with high impacts, such as [58] and [18], which successfully develop Markov chain Monte Carlo (MCMC) and stochastic variational inference (SVI) techniques for system identification. More recently, the rise of variational auto-encoders [59] has enabled a new and powerful view on inference techniques [60] of hybrid systems. A distinct property of such approaches is their reliance on endto-end differentiability and the need to relax discrete variables in order to perform inference.
In the domain of learning-for-control, the notion of switching systems is directly related to the paradigm of model-free hierarchical reinforcement learning (HRL) [61], [62], which combines simple representations to build complex policies. Here it is useful to differentiate between two concepts of hierarchical learning, namely temporal [63], and state abstractions [64]. In their seminal work [65], [66] the authors build on the framework of semi-Markov decision processes (SMDP) [67] to learn activation/termination conditions of temporally extended actions (options) for solving discrete environments. Additionally, pioneering work in optimizing hierarchical control structures with temporally extended actions for robotic applications is developed in [68] specifically and [69] more broadly. Further recent work has focused on different formulations of the option framework that facilitate simultaneous discovery and learning of options [70]- [74].
However, the concept of state abstraction -the aggregation of state-action spaces into sub-regions, each governed by local dynamics and control -carries the most apparent parallels to the classical view of hybrid systems. In [27], a proof of convergence for RL in tabular environments with state abstraction is presented, while [75] does a comprehensive study of different abstraction schemes and gives a formal definition of the problem. Furthermore, recent work has shown promising results in solving complex tasks by combining local linear policies, albeit while still leveraging a complex neural network architecture as an upper-level policy [76].
Switching systems serve as a powerful tool in behavioral cloning. For example, [77] combines hidden Markov models (HMMs) with Gaussian mixture regression to represent trajectory distributions. In contrast, [72] uses a semi-hidden Markov model (HSMM) to learn hierarchical policies, and [78] introduces switching density networks for system identification and behavioral cloning. Finally, a fully Bayesian framework for the hierarchical decomposition of policies is presented in [79], albeit while considering known transition dynamics.
In light of the reviewed literature, we highlight the main differences that distinguish this paper from the approaches mentioned above. First, this work leverages probabilistic hybrid dynamic networks as hierarchical representations of nonlinear open-and closed-loop behaviors. Contrary to standard piecewise autoregressive exogenous systems (PWARX), HDBN can easily integrate stochasticity and nonlinear switching boundaries, leading to more refined and less redundant segmentations. Furthermore, by pursuing an abstraction over the states instead of time, we circumvent the need to infer so-called termination policies, a characteristic of the option framework. Finally, the proposed hybrid policy search approach formulates a non-convex infinite horizon objective that optimizes a hierarchical local polynomial approximation of the value function. This approximation is used to derive stationary switching feedback controllers. In contrast, trajectory optimization and model predictive control techniques are often cast as sequential convex programs and optimize a fixed horizon objective that yields time-variant value functions and controllers.

III. PROBLEM STATEMENT
Consider the discrete-time optimal control problem of a stochastic nonlinear dynamical system to be defined as an infinite horizon Markov decision processes (MDP). An MDP is an abstraction of an environment defined over a state space X ⊆ R d and an action space U ⊆ R m . The probability of a state transition from state x to state x' by applying action u is governed by the Markovian time-independent density function p(x'|x, u). The reward r(x, u) is a function of the state x and action u and is discounted over time by a factor ϑ ∈ [0, 1). The state-dependent policy π(u|x), from which the actions are drawn, is a density determining the probability of an action u given a state x. The general objective in an infinite horizon optimal control problem is to maximize the expected cumulative sum of discounted rewards V π (x) = E ∞ t ϑ t r , where V π denotes as the state-value function under the policy π, starting form an initial state distribution µ 1 (x).
Given the context of this work and our choice to model the system with switching linear models, we introduce to the MDP formulation a new hidden discrete variable z, an indicator of the currently active local regime. The resulting transition dynamics can then expressed by a factorized density Fig. 2: Graphical model of recurrent autoregressive hidden Markov models (rARHMMs) extended to support hybrid controls. rARHMMs are hybrid dynamic Bayesian networks that explicitly allow the discrete state z to depend on the continuous variables x and u, as highlighted in red.
function p(x', z'|x, u, z) = p(z'|z, x, u)p(x'|x, u, z'), which we depict as a graphical model in Figure 2 and discuss in further detail in the upcoming section. In the same spirit of simplification through hierarchical modeling, we employ a mixture of switching polynomial controllers π(u|x, z), associated with local polynomial value functions V π (x, z). The resulting framework becomes a combination of filtering to infer the active local dynamics denoted by z and optimal control to apply appropriate actions u given x and z.

IV. HYBRID DYNAMIC BAYESIAN NETWORKS
We focus on recurrent autoregressive hidden Markov models (rARHMMs) as a representation of closed-loop stochastic hybrid systems. The rARHMM is a special case of recurrent switching linear dynamical systems (rSLDS) [18], also known as augmented SLDS [80]. In contrast to rSLDS, an rARHMM lacks an observation model and directly describes the internal state with an additive noise process. We extend rARHMMs to support exogenous and endogenous inputs in order to simulate the open-and closed-loop behaviors of driven dynamics. Figure 2 depicts the corresponding graphical model, which closely resembles the linear boundary PWARX.
An rARHMM with K regions models the trajectory of a hybrid system as follows. The initial continuous state x 1 ∈ R d and continuous action u 1 ∈ R m are drawn from a pair of Gaussian and conditional Gaussian distributions 1 , respectively. The initial discrete indicator z 1 is a one-hot random vector modeled by a categorical density parameterized by ϕ The transition of the continuous state x t+1 and actions u t are modeled by linear-Gaussian dynamics where (A, B, c, K, Ω, Λ, ∆) are matrices and vectors of appropriate dimensions w.r.t. x and u. φ(x) are polynomial state features of arbitrary degree. The discrete transition probability p(z t+1 |z t , x t , u t ) is governed by K categorical distributions parameterized by state-action dependent multi-class logistic functions where f may have any type of x and u features, e.g., polynomial or neural. A set of vectors ω ij parameterize all transition combinations i → j ∀i, j ∈ [1, K]. Figure 1 depicts realizations of different transition functions that lead to a variety of state space decompositions. This representation of switching dynamics has a significant advantage over other non-recurrent hybrid models [58], [81], since it couples discrete and continuous dynamics of an HMM in both directions. This aspect has significant implications for the model's expressiveness and ability to capture the underlying dynamics of interesting physical applications, as it limits redundancies in the hierarchical decomposition of the state-action space. For example, one may consider a case in which one component can explain the local dynamics in the neighborhood of multiple non-connected discrete states. To achieve a sharp decision boundary in such scenarios, a nonrecurrent model has to duplicate the continuous dynamics with a different transition probability for each set of neighboring regions. This duplication leads to redundant discrete states with the same continuous dynamics while differing in their switching behavior. One way to circumvent this explosion is to consider a hierarchical abstraction over meta regions similar to factorial hidden Markov model (FHMM) [82]. Nonetheless, such representations may still require multiple hierarchy levels to match the expensiveness of a recurrent transition, which compactly parameterizes a continuum of transition probabilities.
The remainder of this paper focuses on using these hybrid models in three ways. First, an open-loop setting that treats the control u as an exogenous input is used for automatically identifying nonlinear systems via decomposition into continuous and discrete switching dynamics. Second, a closedloop setting that assumes the control u to originate from a nonlinear controller. We show that this setting can simultaneously decompose dynamics and control in an behavioral cloning scenario. Finally, we leverage the same framework in a model-based hybrid reinforcement learning algorithm to learn switching controllers of general nonlinear systems.

V. BAYESIAN INFERENCE OF HYBRID SYSTEM IDENTIFICATION AND BEHAVIOR CLONING
In this section, we sketch the outline of an expectationmaximization/Baum-Welch algorithm [84]- [86] for inferring the parameters θ of an rARHMM given time-series observations (X, U). The resulting algorithm can be used twofold. First, it can be applied to perform automatic hybrid system identification to learn the open-loop dynamics of nonlinear systems given state-action observations. Second, it can clone the closed-loop behavior of a nonlinear controller and decompose it into a set of local experts.
Our developed approach is related in some aspects to the Baum-Welch algorithms proposed in [87] and [72]. However, we introduce suitable priors over all parameters and derive a maximum a posteriori (MAP) technique with a stochastic maximization step and hyperparameter optimization. In our experience, the priors and noisy gradient estimate significantly reduce the sensitivity of EM w.r.t. initialization and appear to be less prone to get stuck in bad local minima -an effect well-studied in neural networks [88].
Although our procedure is not fully Bayesian like methods that rely on Gibbs sampling [18], [58], it has computational and predictive advantages. On the one hand, Gibbs sampling can suffer from slow convergence in high dimensional spaces leading to an overall high computational cost [89]. On the other hand, standard Gibbs sampling-based approaches are not flexible enough to incorporate neural transition predictor functions due to their reliance on conditionally conjugate computation. Finally, a good prior specification is crucial in small data regimes since a vague prior may dominate the predictive posterior and effectively cause under-fitting. We implement a hyperparameter optimization scheme that elevates this concern by optimizing the prior parameters via empirical Bayes [90], thus attenuating the prior influence and improving the predictive performance significantly.

A. Maximum A Posteriori Optimization
Consider again the rARHMM in Figure 2 where the continuous state x and action u are observed variables, while the K-region indicators z are hidden. To infer the model parameters, we assume a dataset consisting of N state-action trajectories D = {(X n , U n )} N n=1 , each of length T , where (X n , U n , Z n ) represent an entire trajectory. The inference objective for system identification is to maximize the log-posterior probability of the observations (X n , U n ) conditioned on the free parameter set θ where p(X n , U n , Z n |θ) is the complete-data likelihood of a single trajectory and factorizes according to and p(θ|h) is the factorized parameter prior We choose all priors to be conjugate or semi-conjugate w.r.t. their likelihoods. Therefore, we place a normal-Wishart (NW) prior on the initial state distribution (µ k , Ω k ) ∼ NW(0, κ 0 , Ψ 0 , ν 0 ), and a matrix-normal-Wishart (MNW) on the linear transition dynamics . The initial discrete state takes a Dirichlet prior ϕ ∼ Dir(τ 0 ), while the logistic transition parameters are governed by a non-conjugate zero-mean Gaussian prior with diagonal precision ω ik ∼ N(0, αI). Finally, we place a separate matrix-normal-Wishart prior on the conditional action likelihood The choice of priors is not restricted to these distributions. Depending on modeling assumptions, one can assume dynamics with diagonal noise matrices and pair them with gamma distribution priors. Moreover, if the system is known to have a state-independent noise process, the K Wishart and gamma priors can be tied across components, leading to a more structured representation.

B. Baum-Welch Expectation-Maximization
On closer examination of Equations (1) and (2), we observe that the optimization problem is non-convex with multiple local optima, since the complete-data likelihood p(D, Z|θ) = p(X, U, Z|θ) can follow complex multi-modal densities. Another technical difficulty is the summation over all possible trajectories of the hidden variables Z n , which is of computational complexity O(N K T ) and is intractable in most cases. Expectation-maximization algorithms overcome the latter problem by introducing a variational posterior distribution over the hidden variables q(Z n ) and deriving a lower bound on the complete log-probability function We can find a point estimate of the parameters θ MAP by following a modified scheme of EM, alternating between an expectation step (E-step), in which the lower bound in Equation (3) is maximized w.r.t. the variational distributions q(Z n ) given a parameter estimateθ, a maximization step (M-step), that updates θ given (q(Z n ),ĥ), and finally an empirical Bayes step (EB-step) that updates h given (q(Z n ),θ). A sketch of the overall iterative procedure is presented in Algorithm 1.

1) Exact Expectation
Step: Maximizing the lower bound w.r.t. q(Z n ) is determined by reformulating Equation ( This form of the lower bound implies that the optimal variational distributionq(Z n ) minimizes the Kullback-Leibler divergence (KL), meaninĝ Algorithm 1: Expectation-Maximization for System Identification and Behavioral Cloning input: K, X, U, h This update tightens the bound if the posterior modelq(Z n ) belongs to the same family of the true posterior [16]. Notice that the E-step is independent of the prior p(θ). Moreover, Equation (4) indicates that the E-step reduces to the computation of the smoothed marginals p(z n t |x n 1:T , u n 1:T ,θ) under the current parameter estimateθ. Following [84] and [91], we derive a forward-backward algorithm, which enables closedform and exact inference: γ n 2) Stochastic Maximization Step: After performing the Estep and computing the smoothed posteriors, we are able to evaluate the lower bound and maximize it w.r.t. θ given (q(Z n ),ĥ). By plugging Equation (4) into (3) and leveraging conditional independence, we arrive at the complete logprobability function The optimization of Q is commonly done via coordinate ascent. Simpler models, e.g., Gaussian-and Binomial-HMMs, lead to an exact convex M-step with closed-form optimality conditions. This is not the case in rARHMM, given the choice of nonlinear transition predictor functions. Such a choice leads to an approximate, non-convex M-step that requires gradient-based updates. In this case, stochastic optimization is recommended [92] as batched noisy gradient estimates allow algorithms to escape shallow local minima and reduce the computational cost that comes with evaluating the gradients for all data instances. Consequently, when implementing the M-step, we apply stochastic optimization on the transition parameters ω. We use a stochastic gradient ascent direction with an adaptive learning rate ε and batch size M [92] For the parameters with conjugate priors, we derive closedform optimality conditions. Effectively, this part of the optimization constitutes formulating the posterior distribution of all remaining parameters and taking the mode of each posterior density for a point estimate update.
By considering only relevant terms, we write the optimization of the initial gating parameter ϕ as Analogously, the objective terms related to the dynamics parameter (A k , B k , c k , Λ k ) are also decoupled to k parts and, finally, to learn closed-loop behavior, we infer the con- Due to space constraints, we will refrain from stating the explicit solution for these former problems. Instead, we provide a general outline of how to compute these posteriors and their modes in Appendix I and II.
3) Approximate Empirical Bayes: Inference techniques that leverage data-independent assumptions run the risk of prior miss-specification. In our MAP approach, the priors are weakly informative and carry little information. Their main purpose is to regularize greedy updates that might lead to premature convergence. However, when there is little data, the priors, especially those on the precision matrices 3 , may dominate the posterior probability, leading to over-regularization and under-fitting of the objective. Empirical Bayes approaches remedy this issue by integrating out the parameters θ and optimizing the marginal likelihood w.r.t. the hyperparameters h [90]. In our setting, marginalizing all hidden quantities does not admit a closed-form formula. An approximate approach to empirical Bayes is to interleave the E-and M-steps with hyperparameter updates that optimize the lower bound given an estimate of parametersθ and a step size The necessary gradients for performing the empirical Bayes updates are provided in Appendix II.

VI. HIERARCHICAL REINFORCEMENT LEARNING
FOR HYBRID DYNAMICAL SYSTEMS In our review of related work in Section II, we highlight that many successful classical hybrid control algorithms are often limited to computationally expensive trajectory-centric control policies, e.g. model predictive control. This type of controller is disadvantageous in applications that require fast reactive feedback control with broad coverage over the state-action space. Another common drawback of classical hybrid control approaches is their reliance on linear separation boundaries of the local dynamics. On the other hand, we mentioned several RL-based approaches that learn global time-invariant policies. Nonetheless, these algorithms are exclusively model-free and mostly rely on time abstraction and temporally extended actions, which requires learning task-specific termination policies, in addition to the local controllers.
In this section, we present a stochastic infinite horizon optimization technique that leverages the structure and properties of hybrid systems under the paradigm of state abstraction. Our algorithm extends the step-based formulation of relative entropy policy search (REPS) [94]- [96] by introducing a discrete state variable z and taking into account the structure of hybrid dynamics. Our approach, hybrid REPS (Hb-REPS), leverages the state-action-dependent nonlinear switches p(z'|z, x, u) as a task-independent upper-level coordinator to a mixture of K lower-level stationary polynomial policies π(u|x, z). While the proposed approach shares many features with [72], our formulation relies on a state-abstraction representation of hybrid systems and embeds the hierarchical model structure into the optimization problem in order to learn a hierarchy over the global value function. In contrast, [72] operates in the framework of semi-Markov decision processes and optimizes a mixture over termination and feedback policies without considering the existence of a hierarchical structure in the space of dynamics and value functions.

A. Infinite-Horizon Stochastic Hybrid Control
In the REPS framework an optimal control problem is presented as an iterative trust-region optimization for a discounted average-reward objective under a stationary stateaction distribution π(u|x, z)µ(x, z), Equation (5a). The trustregion is formulated as a KL [97], Equation (5c). Its purpose is to bound the information loss between iterations. The REPS formulation explicitly incorporates a dynamics consistency constraint, Equation (5b), that describes how the stochastic state of the system evolves. This property motivated our algorithmic choice when integrating the hybrid dynamics. The following describes the optimization solved during a single iteration of what we refer to as hybrid REPS where µ(x, z) is the stationary mixture distribution, q(x, u, z) is the trust-region reference distribution, and the constraint in Equation (5d) guarantees the normalization of the state-action distribution. The factor 1 − ϑ represents the probability of an infinite process to reset to an initial distribution µ 1 (x, z). The notion of resetting is necessary to ensure ergodicity and allows the interpretation of ϑ as a discount factor and regularization of the MDP [96], [98].

B. Optimality Conditions and Dual Optimization
For convenience, let p(x, u, z) = π(u|x, z)µ(x, z) for which the marginalization identity p(x, u, z) du = µ(x, z) holds. Using the method of Lagrangian multipliers [99], we can solve the constrained primal problem w.r.t. state-action distribution p(x, u, z) where η is the Lagrangian variable associated with Equation (5c), and A(x, u, z) the advantage function given as The function V (x, z) is the state-value function, which appears naturally in REPS as the Lagrangian function associated with Equation (5b). The full Lagrangian of the primal can be found in Appendix III. By substituting the optimal parameter p * back into the Lagrangian and factorizing q(x, u, z), we arrive at the dual function G as a function of the remaining Lagrangian variables η and V where q(z|x, u) is the posterior over z given the observations x and u. In Section V, we derived a forward-backward algorithm for inferring these probabilities, allowing us to compute the expectation over z. The expectations over X and U are usually analytically intractable, therefore, we approximate them given samples from the reference distribution q(x, u). The multipliers η and V are then obtained by numerically minimizing the dual G(η, V ) arg min that acts as the upper bound on the primal objective.

C. Stationarity of State Distribution Mixtures
The dynamics Equation (5b), which ensures the stationarity of µ(x, z), uncovers an interesting aspect of our optimization problem. A careful inspection of that integral equation reveals that a mixture distribution µ(x, z) propagated through the mixture dynamics p(x', z'|x, u, z) results in a mixture µ(x', z') that keeps growing with every pass, leading to an explosion in the number of components of the joint state distribution. This problem highlights a crucial computational issue of trajectorybased optimal control approaches of stochastic hybrid systems, as it becomes expensive to maintain the full mixture state distribution µ(x, z) after several time steps. Common solutions to this issue usually involve falling back to crude Gaussian mixture reduction techniques [100] that inadvertently sacrifice information and blur the distribution as it progresses in time.
This consideration has motivated us to formulate the optimal control problem as one that focuses on finding a stationary solution for µ(x, z). We hypothesize, that including Equation (5b) as a constraint leads to an optimal mixture policy π(u|x, z) that acts as a filter and dampens certain modes of µ(x', z') with little contribution to the average-reward objective. However, we recognize that a more in-depth analysis of the nature of this integral equation is necessary.

D. Modeling Dynamics and State-Value Function
Up to this point, our derivation has been generic. We have made no assumptions on initial distributions µ 1 (x, z), the dynamics p(x', z'|x, u, z), or the value function V (x, z). Now we introduce the local Gaussian linear dynamics and logistic switching described in Section IV and assume these representations to be available in parametric form as a result of separate learning process. Furthermore, we model the statevalue function with local n-th degree polynomial functions is the state-feature vector which contains all polynomial features of the state x, and ω z is the parameter vector assigned to the different regions.
Under these assumptions, we can leverage the available joint density µ 1 (x, z) and p(x'|x, u, z) to compute the necessary expectations in Equation (7) This computation allows our approach to capture the stochasticity of the dynamics and delivers an estimate of the advantage function A(x, u, z) instead of the temporal difference (TD) error in the general REPS framework [94]. Ultimately, this leads to better estimates of the expected discounted future returns.
Practically, these integrals can be either naively approximated by applying Monte Carlo integration [101], or, more efficiently, by leveraging the structure of the integrand V (x', z'), and using Gauss-Hermite cubature rules for exact integration over polynomial functions [102].

E. Maximum-A-Posteriori Policy Improvement
Another important advantage of this model-based approach becomes evident when considering the policy improvement step in the REPS framework. The policy update is incorporated into the optimality condition of the stationary state-action distribution p(x, u, z) = π(u|x, z)µ(x, z) in Equation (6). As a consequence, updating the mixture policies π(u|x, z) requires the computation of state probabilities µ(x, z), which in turn require knowledge of the dynamics model. This issue is circumvented in other model-free realizations of REPS by introducing a crude approximation to enable a model-free policy update nonetheless. In [103], the authors postulate that the new state distribution µ(x, z) is usually "close enough" to the old distribution q(x, z), thus allowing the ratio q(x, z)/µ(x, z) to be ignored when a weighted maximum-likelihood fit of the actions u is performed to update π. initialize: q(u|x, z), ω z , η While the assumption of "closeness" may be practical and empowers many successful flavors of REPS, it is crucial to be aware of its technical ramifications, as it undermines the primary motivation of a relative entropy bound on the stateaction distribution in Equation (5c). This aspect is unique in the REPS framework when compared to other state-of-the-art approximate policy iteration algorithms [7], [9], [104], that optimize a similar objective, albeit with a relaxed bound that only limits the change of the action distribution π.
In contrast, our proposed algorithm leverages the modeled continuous-discrete dynamics and updates the policy π(u|x, z) with the correct weighting. The optimality condition in Equation (6) is satisfied by deriving a weighted maximum a posteriori estimate based on the state-action distribution p(x, z, u), and implicitly updating π(u|x, z). This procedure is equivalent to a modified Baum-Welch expectationmaximization algorithm that fits a closed-loop rARHMM, as derived in Section V. The difference is that the EM objective in Equation (1) has to be augmented with the importance weights from Equation (6) arg max θ log N n=1 z n w n p(X n , U n , Z n |θ)p(θ), This augmentation leads to weighted M-and EB-steps while the E-step is not altered.
Note that during the policy improvement step, we can either assume an a priori estimate of the open-loop dynamics p 0 (x', z'|x, u, z) and only update the control parameters corresponding to the conditional π(u|x, z), or we can continuously update p(x', z'|x, u, z) as more data iteratively becomes available. A compact sketch of the overall optimization process is available in Algorithm 2.

VII. EMPIRICAL EVALUATION
In this section, we benchmark different aspects related to the inference of rARHMMs and the learning of hybrid controllers    Figure 3.
The values reflect the total number of parameters of each model. The values in parentheses represent the hidden layer sizes S of the neural models and the number of discrete components K for the (r)ARHMM, respectively.
via Hb-REPS. First, we assess the predictive performance of rARHMMs in an open-loop setting to validate the choice of this representation of hybrid dynamics. Second, we test the inference on closed-loop rARHMMs and their ability to capture and decompose an expert nonlinear controller in a behavioral cloning scenario. Finally, we use rARHMMs in the proposed RL algorithm Hb-REPS to solve the infinite horizon stochastic control objective and optimize piecewise local polynomial controllers and value functions. Here we make no claim to the absolute efficiency of our approach when compared to other state-of-the-art algorithms. Instead, we aim to provide an empirical proof-of-concept that supports further research into sample-based hybrid system optimization as a framework for structured nonlinear control.

A. Hybrid System Identification Examples
We empirically benchmark rARHMMs on nonlinear systems. We aim to quantify the quality of learned open-loop models and their ability to capture the underlying dynamics. For this purpose, we set up a direct comparison to popular representations for dynamics in a long-horizon and limited-data setting.Public code base available on https://github. com/hanyas/sds.
In this evaluation, we focus on rARHMMs with exogenous inputs. We learn the dynamics of three simulated deterministic systems; a bouncing ball, an actuation-constrained pendulum, and a cart-pole system. We compare the predictive accuracy of rARHMMs to classical non-recurrent autoregressive hidden Markov models (ARHMMs) [17], feed-forward neural nets (FNNs), Gaussian processs (GPs) 4 , long-short-term memory networks (LSTMs) [105], and recurrent neural networks (RNNs). During the evaluation, we collected segregated training and test datasets. The training dataset is randomly split into 24 groups, each containing a subset of trajectories, and used to train different instances of all representations. These instances are then tested on the test dataset. All neural models have 2 hidden layers, which we test for a variety of different layer sizes, S = {16, 32, 64, 128, 256, 512} for FNNs, S = {16, 32, 64, 128, 256} for RNNs, and S = {16, 32, 64, 128} for LSTMs. In the case of (r)ARHMMs, we test for different number of component K, dependent on the task. As a metric, we evaluate the NMSE, averaged over the 24 data splits for a range of horizons. During evaluation, we comb through the test trajectories step by step and predict the given horizon. Moreover, in Table I,  representations in terms of their total number of parameters. a) Bouncing Ball: This example is a canonical instance of a dual-regime hybrid system due to the hard velocity switch at the moment of impact. We simulate the dynamics with a frequency of 20 Hz and collect 25 training trajectories with different initial heights and velocities, each 30 s long. This dataset is split 24 folds with ten trajectories, 10 × 150 data points, in each subset. The test dataset consists of 5 trajectories, each 30 s long. We evaluate the NMSE for horizons h = {1, 20, 40, 60, 80} time steps. We did not evaluate a GP model in this setting due to the long prediction horizons that led to a very high computational burden. The (r)ARHMMs are tested for K = 2. The logistic link function of an rARHMM is parameterized by a neural net with one hidden layer containing 16 neurons. The results in Figure 3 show that the rARHMM approximates the dynamics well and outperforms both ARHMMs and the neural models. The logistic link function of the rARHMM is parameterized by a neural net with one hidden layer containing 24 neurons. As shown in Figure 3, the forecast evaluation provides empirical evidence for the representation power of rARHMMs in both smooth and discontinuous state spaces. FNNs and GPs perform well in the smooth Cartesian observation space, and struggle in the discontinuous space, similar to RNNs and LSTMs.

B. Hierarchical Closed-Loop Behavioral Cloning
Before applying the RL method proposed in Section VI-A, we first want to analyze the closed-loop rARHMM with endogenous inputs as a behavioral cloning framework. The task is to reproduce the closed-loop behavior of expert policies on challenging nonlinear systems. For this purpose, we train two different feedback experts on the pendulum and cartpole. The two environments are simulated at 50 Hz and are influenced by static Gaussian noise with a standard deviation σ = 1 × 10 −2 . The experts are two-layer neural nets with 4545 parameters (pendulum) and 17537 parameters (cartpole), optimized with the SAC algorithm [9]. For cloning, we construct two 5-regime rARHMMs with local polynomial policies of the third order. The hybrid controllers have a total number of parameters of 100 (pendulum) and 280 (cart-pole). Learning is realized on a dataset of 25 trajectories, each 5 s long, for each environment and using the EM technique from Section V. The decomposed controllers complete the task of swinging up and stabilizing both systems with over 95% success rate. Figure 4 shows the phase portraits of the unforced dynamics and closed-loop control identified during cloning. Figure 5 depicts sampled trajectories of the hybrid policies highlighting the switching behavior.

C. Reinforcement Learning for Hybrid Systems
Finally, we evaluate the qualitative performance of the proposed hybrid policy search algorithm Hb-REPS on two nonlinear stochastic dynamical systems: an underpowered pendulum swing-up and a cart-pole stabilization task that explicitly simulates an abrupt switch in dynamics when the cart hits an elastic wall.
We compare the performance of Hb-REPS to two baselines. The first is a vanilla version of REPS that does not maintain any hierarchical structure and uses nonlinear function approximators with random Fourier features (RFFs) [106] to represent both policy and value function. The second baseline assumes a hierarchical policy structure and a nonlinear value function with Fourier features. This baseline is akin to what is implemented in [72], albeit with a hierarchy based on state abstraction rather than time. We will refer to this algorithm version as hierarchical REPS (Hy-REPS). We assume an offline learning phase in which the hybrid models are learned from pre-collected data. a) Pendulum Swing-up: In this experiment, the powerlimited pendulum is simulated at 50 Hz and perturbed by Gaussian noise with a standard deviation σ = 1 × 10 −2 . The REPS agent relies on a policy and value function with 50 and 75 Fourier features, respectively. Hy-REPS assumes a similar form of the value function but with five third-order polynomial policies. Hb-REPS represents both policy and value function with five third-order polynomials. Empirical results in Figure 7 (left half) feature comparable learning performance of all algorithms over ten random seeds. Every iteration involves 5000 interactions with the environment. We provide a phase portrait of the closed-loop behavior for a qualitative assessment of the final stationary hybrid policy. b) Cart-pole Stabilization: This evaluation features a cartpole constrained by an elastic wall modeled as a spring. The dynamics are linearized around the upright position, naturally resulting in a two-regime hybrid system. The environment is simulated at 100 Hz and perturbed by Gaussian noise with a standard deviation σ = 1 × 10 −4 . The REPS policy and value function both use 25 random Fourier features. Hy-REPS adopts the same value function structure with two affine linear policies. Hb-REPS also assumes two affine linear policies combined with two second-order local value functions. The results in Figure 7 (right half) depict matching learning performance of the three approaches over ten random seeds. Every iteration involves 2500 interactions with the environment.

VIII. DISCUSSION
We presented a data-driven view of hybrid system identification and control that serves as an alternative paradigm to common popular techniques. Our approach is not restricted to the class of explicit hybrid dynamics and can be seen as a general approach to structured identification and nonlinear control. We argue that this structure often exists under the hood of complex neural representation. Therefore, making it explicit may offer an avenue to apply Occam's razor and regularize over-parameterized representations. Initial empirical results support this motivation. The proposed hybrid reinforcement learning technique can do without neural networks and instead rely on a hierarchy of local polynomial representations dictated by a hierarchical structure.
Nonetheless, the application of this work is limited to simple low-dimensional systems. Although a viable alternative to expensive mixed-integer optimization, the inference techniques used in this paper still present a bottleneck in the face of scalability to more complex systems and higher dimensions. While our MAP approach significantly improves the quality of expectation-maximization solutions, it nevertheless struggles in more challenging environments.
A possible course of action is to investigate Bayesian nonparametric extensions of hybrid dynamic Bayesian networks based on non-conjugate variational inference. Fully Bayesian methods tend to improve learning in large structured models significantly. Another potential avenue of research is to improve the hybrid reinforcement learning framework by considering the control-as-inference paradigm. Such approaches may offer ways of integrating the Bayesian structure of the models into the control optimization and, in the process, achieve an uncertainty-aware approach that is better equipped to deal with the exploration-exploitation dilemma.

APPENDIX I EXPONENTIAL FAMILY
Our work mainly considers random variables with probability density functions belonging to the exponential family. The unified minimal parameterization of this class of distributions lends itself for convenient and efficient posterior computation when paired with conjugate priors.
We assume the natural form for a probability density of a random variable x where h(x) is the base measure, η are the natural parameters, t(x) are the sufficient statistics and a(η) is the log-partition function, or log-normalizer. Following the same notation, a conjugate prior g(η|λ) to the likelihood f (x|η) has the form with prior sufficient statistics t(η) = [η, −a(η)] and hyperparameters λ = [α, β] . By applying Bayes' rule, we can directly infer the posterior q(η|x) where the posterior natural parameters ρ(x, λ) are a function of the likelihood sufficient statistics t(x) and prior hyperparameters [α, β] ρ(x, λ) = α + t(x), β + 1 .
The structure of the resulting posterior reveals a simple recipe for data-driven inference. By moving into the natural space, the posterior parameters are computed by combining the prior hyperparameters with the likelihood sufficient statistics and log-partition function. By definition, every exponential family distribution has a minimal natural parameterization that leads to a unique decomposition of these quantities [107].

APPENDIX II CONJUGATE POSTERIORS
We present an outline of all M-step updates. We use an adapted form of the exponential natural parameterization, as it offers a clear methodology on how to derive and implement such updates for all relevant distributions.

A. Categorical with Dirichlet Prior
A weighted categorical likelihood over a one-hot random variable z with size K has the form . . .   Finally, the gradients of prior parameters needed during the empirical Bayes incremental updates are ∂ log p(ϕ) ∂τ k = ψ( K k=1 τ k ) − ψ(τ k ) + log ϕ k , where ψ stands for the digamma function.

C. Linear-Gaussian with Matrix-Normal-Wishart Prior
A weighted linear-Gaussian likelihood takes a random variable x ∈ R d and returns a random variable y ∈ R m according to a linear mapping A : R d → R m p(Y|X, A, V) = N n=1 N(y n |x n , A, V) wn ∝ exp where w n are the weights and W = diag(w n ) is the diagonal weight matrix. The data matrices X and Y are of size d × N and m × N respectively. The conjugate prior p(A, V) is a matrix-normal-Wishart with zero mean The posterior q(µ, Λ) is matrix-normal-Wishart q(A, V) = N(A|M, V, K) W(V|Ψ, ν) The mode mapping and precision of a matrix-normal-Wishart areÂ = M andΛ = (ν − m)Ψ, respectively. The standard posterior parameters are The empirical Bayes gradients of the log-prior are