Non-smooth Dynamical Systems Research Articles

Solving Mathematical Programs with Complementarity Constraints Arising in Nonsmooth Optimal Control

AbstractThis paper examines solution methods for mathematical programs with complementarity constraints (MPCC) obtained from the time-discretization of optimal control problems (OCPs) subject to nonsmooth dynamical systems. The MPCC theory and stationarity concepts are reviewed and summarized. The focus is on relaxation-based methods for MPCCs, which solve a (finite) sequence of more regular nonlinear programs (NLP), where a regularization/homotopy parameter is driven to zero. Such methods perform reasonably well on currently available benchmarks. However, these results do not always generalize to MPCCs obtained from nonsmooth OCPs. To provide a more complete picture, this paper introduces a novel benchmark collection of such problems, which we call . The problem set includes 603 different MPCCs and we split it into a few representative subsets to accelerate the testing. We compare different relaxation-based methods, NLP solvers, homotopy parameter update and relaxation parameter steering strategies. Moreover, we check whether the obtained stationary points allow first-order descent directions, which may be the case for some of the weaker MPCC stationarity concepts. In the best case, the Scholtes’ relaxation (SIAM J. Optim. 11, 918–936, 2001) with (Math. Program. 106, 25–57, 2006) as NLP solver manages to solve 73.8% of the problems. This highlights the need for further improvements in algorithms and software for MPCCs.

Fractional Differential Equations with Impulsive Effects

This paper discusses impulsive effects on fractional differential equations. Two approaches are taken to obtain our results: either with fixed or changing lower limits in Caputo fractional derivatives. First, we derive an existence result for periodic solutions of fractional differential equations with periodically changing lower limits. Then, the impulsive effects are modeled for fractional differential equations regarding the nonlinearities rather than the initial value conditions. The proposed impulsive model differs from common discontinuous and nonsmooth dynamical systems.

Well-posedness and optimal control of a nonsmooth fractional dynamical system

Well-posedness and optimal control of a nonsmooth fractional dynamical system

Finite Elements with Switch Detection for numerical optimal control of nonsmooth dynamical systems with set-valued heaviside step functions

This paper develops high-accuracy methods for the numerical solution of optimal control problems subject to nonsmooth differential equations with set-valued Heaviside step functions. An important subclass of these systems are Filippov systems. By writing the Heaviside step function as the solution map of a linear program and using its optimality conditions, the initial nonsmooth system is rewritten into an equivalent Dynamic Complementarity System (DCS). The Finite Elements with Switch Detection (FESD) method (Nurkanović et al., 2024) was originally developed for Filippov systems transformed via Stewart’s reformulation into DCS (Stewart, 1990). This paper extends it to the above mentioned class of nonsmooth systems. The key ideas are to start with a standard Runge–Kutta method for the DCS and to let the integration step sizes to be degrees of freedom. Then, additional conditions are introduced to allow implicit but accurate switch detection and to remove possible spurious degrees of freedom when no switches occur. The theoretical properties of the FESD method are studied. The motivation for these developments is to obtain a computationally tractable formulation of nonsmooth optimal control problems. Numerical simulations and optimal control examples are used to illustrate the favorable properties of the proposed approach. All methods introduced in this paper are implemented in the open-source software package nosnoc (Nurkanović and Diehl, 2022).

Computing subgradients of convex relaxations for solutions of parametric ordinary differential equations

A novel subgradient evaluation method is proposed for nonsmooth convex relaxations of parametric solutions of ordinary differential equations (ODEs) arising in global dynamic optimization, assuming that the relaxations always lie strictly within interval bounds during integration. We argue that this assumption is reasonable in practice. These subgradients are computed as the unique solution of an auxiliary parametric affine ODE, analogous to classical forward/tangent sensitivity evaluation methods for smooth dynamic systems. Unlike established subgradient evaluation approaches for nonsmooth dynamic systems, this new method does not require smoothness or transversality assumptions, and is compatible with existing subgradient evaluation methods for closed-form convex functions, as implemented in subgradient evaluation software such as EAGO.jl and MC ++. Moreover, we show that a subgradient for a lower-bounding problem in global dynamic optimization can be directly evaluated using reverse/adjoint sensitivity analysis, which may reduce the overall computational effort for an overarching global optimization method. Numerical examples are presented, based on a proof-of-concept implementation in Julia.

Nonlinear characteristics identification of an impact oscillator with a one-sided elastic constraint

Impacting systems are widely used in many engineering applications, such as self-propelled robots, energy harvesting and percussive drilling, which exhibit rich and complex nonlinear phenomena. Among these applications, predicting nonlinearities and estimating system parameters are of great interest of nonlinear dynamics research community. Backbone curve is an analytical tool that captures the frequency–amplitude dependence of nonlinear systems. In this paper, we estimate the impacting stiffness of a single-degree-of-freedom non-smooth dynamical system qualitatively by using the backbone curve. It was found that an increase of the impacting stiffness may lead to lowering the backbone curve. An adaptive differential evolution algorithm with the Metropolis criterion is proposed to identify the parameters of the impacting system quantitatively using experimental data, which are consistent with our theoretical predictions. Finally, the identified parameters are verified, and the limitations of the backbone curve are drawn. The nonlinear characteristics identification method studied in this paper can be extended to other vibro-impact systems and is potentially applicable to structural health monitoring and robotic sensing.

A nonlinear transient-dynamics approach to atopic dermatitis: Role of spontaneous remission

Atopic dermatitis (AD) is a common skin disease that can occur in all age groups. An intriguing phenomenon associated with AD is spontaneous remission, in which the symptoms can disappear even without any treatment, especially for patients at a young age. From the point of view of dynamical evolution, spontaneous remission in AD is a transient phenomenon. A clinic implication is that, if the transient time is short, then aggressive treatment may not be necessary. A key question is thus, statistically, how long the transient time can be? Due to the lack of clinic data, mathematical modeling is a viable approach to addressing this question. Modeling AD as a nonsmooth dynamical system and regarding the disease as a transient phenomenon with spontaneous remission marking the end of the transient, we obtain a quantitative understanding of the statistical characteristics of AD. In particular, we find that, depending on the immune state, two different types of transient behaviors can arise. For type-I spontaneous remission characterized by a healthy immune level with its skin state exhibiting mild oscillations, the transient time is short, which typically occurs in patients between infancy and childhood. In contrast, type-II spontaneous remission is characterized by a low immune level and its skin state exhibits severe oscillations with a long recovering time. Quantitatively, a scaling relation exists between the average transient time (or recovery time) and some key physiological parameters, revealing that the transient is superpersistent in the sense that its average lifetime can diverge in a drastic way: e∞ as a bifurcation parameter approaches a critical value. In this case, the disease is essentially permanent, thereby requiring and justifying active treatment.

Second-order cone programming for frictional contact mechanics using interior point algorithm

We report experiments of an implementation of a primal–dual interior point algorithm for solving mechanical models of one-sided contact problems with Coulomb friction. The objective is to recover an optimal solution with high accuracy and as quickly as possible. These developments are part of the design of Siconos, an open-source software for modelling and simulating non-smooth dynamical systems. Currently, Siconos uses mainly first-order methods for the numerical solution of these systems. These methods are very robust, but suffer from a linear rate of convergence and are therefore too slow to find accurate solutions in a reasonable time. Our main objective is to apply second-order methods to speed up convergence. In this paper, we focus on solving a relaxed formulation of the initial mechanical model, corresponding to the optimality conditions of a convex quadratic minimization problem with second-order cone constraints. We will present in detail a primal–dual interior point algorithm for solving this type of problem. The main difficulty in implementing this algorithm arises from the fact that at each iteration of the algorithm, a change of variable, called a scaling, has to be performed in order to guarantee the non-singularity of the linear system to be solved, as well as to recover a symmetric system. Although this scaling strategy is very nice from a theoretical point of view, it leads to an enormous deterioration of the conditioning of the linear system when approaching the optimal solution and therefore to all the numerical difficulties that result from it. We will describe in detail the numerical algebra that we have developed in our implementation, in order to overcome these problems of numerical instability. We will also present the solution of the models resulting from the problems of rolling friction, for which the cone of constraints is no longer self-dual like the Lorentz cone.

Global dynamics for a nonsmooth dynamical system induced by wing rock

Global dynamics for a nonsmooth dynamical system induced by wing rock

Grazing-Induced Dynamics of the Piecewise-Linear Chua’s Circuit

In this paper, we consider the piecewise-linear Chua’s circuit, which is well known for its rich variety of bifurcation, chaotic and other nonlinear phenomena. Suitable switching boundaries are introduced based on the piecewise-linear representation of Chua’s diode. In this way, we derive analytical conditions for a grazing bifurcation to occur, when one or two families of periodic orbits have a zero-velocity contact with the switching boundaries. In connection to this phenomenon, we also study the focus-center-limit cycle bifurcation and its implications regarding the system dynamics, from both analytical and numerical points of view. Furthermore, a detailed parametric study of Chua’s circuit is carried out via path-following techniques for nonsmooth dynamical systems, implemented via the continuation software COCO. This study reveals the presence of codimension-one bifurcations of limit cycles, such as those mentioned above, as well as classical (fold and period-doubling) bifurcations. The analysis confirms the presence of coexisting attractors, which are produced by a hysteresis loop induced by the interaction of a fold and a focus-center-limit cycle bifurcation.

Geometric slow–fast analysis of a hybrid pituitary cell model with stochastic ion channel dynamics

To obtain explicit understanding of the behavior of dynamical systems, geometrical methods and slow–fast analysis have proved to be highly useful. Such methods are standard for smooth dynamical systems and increasingly used for continuous, non-smooth dynamical systems. However, they are much less used for random dynamical systems, in particular for hybrid models with discrete, random dynamics. Here we propose a geometrical method that works directly with the hybrid system. We illustrate our approach through an application to a hybrid pituitary cell model in which the stochastic dynamics of very few active large-conductance potassium (BK) channels is coupled to a deterministic model of the other ion channels and calcium dynamics. To employ our geometric approach, we exploit the slow–fast structure of the model. The random fast subsystem is analyzed by considering discrete phase planes, corresponding to the discrete number of open BK channels, and stochastic events correspond to jumps between these planes. The evolution within each plane can be understood from nullclines and limit cycles, and the overall dynamics, e.g., whether the model produces a spike or a burst, is determined by the location at which the system jumps from one plane to another. Our approach is generally applicable to other scenarios to study discrete random dynamical systems defined by hybrid stochastic–deterministic models.

Haar wavelet for computing periodic responses of impact oscillators

Haar wavelets have an attractive property for non-smooth dynamical systems as they are capable of modelling sudden changes because of their local multi-resolution characteristics. In this paper, we applied the Haar wavelet collocation method embedding the segment technique to compute and detect periodic responses of an elastic impact oscillator. Comparisons between dynamical responses computed by a direct numerical simulation using a high accuracy Runge–Kutta algorithm and the proposed method are encouraging. Some key parameters of the impact oscillator were changed to demonstrate the effectiveness of the method. Moreover, the time–frequency index feature for the Haar wavelet coefficients describing impacts was proven to be an additional advantage of this method.

Dynamic tipping and cyclic folds, in a one-dimensional non-smooth dynamical system linked to climate models

We study the behaviour at tipping points close to (smoothed) non-smooth fold bifurcations in one-dimensional oscillatory forced systems with slow parameter drift. The focus is the non-smooth fold in the Stommel 2-box, and related climate models, which are piecewise-smooth continuous dynamical systems, modelling thermohaline circulation. These exhibit non-smooth fold bifurcations which arise when a saddle-point and a focus meet at a border collision bifurcation. By using techniques from the theory of non-smooth dynamical systems we are able to provide precise estimates for the general tipping behaviour close to the non-smooth fold as parameters vary. These are different from the usual tipping point estimates at a saddle–node bifurcation, with advanced tipping apparent in the non-smooth case when compared to the behaviour near smoothed approximations to this fold. We also see rapid, and non-monotone, transitions in the tipping points for oscillatory forced systems close to both non-smooth folds and saddle–node bifurcations due to the effects of both phase changes and non-smoothness. These variations can have implications for the prediction of tipping in climate systems, particularly in close proximity to a non-smooth fold.

Analyzing the backbone curve of piecewise-linear non-smooth systems using a generalized bilinear frequency approximation method

Estimating the resonant frequency of non-smooth systems with contacting elements is critical for the design, control, and health monitoring of various mechanical and aerospace structures. The nonlinearity caused by the intermittent contact excludes the use of efficient linear methods; hence, predicting vibration characteristics of these systems becomes computationally challenging. Furthermore, several numerical and analytical approaches that have been developed for analyzing smooth nonlinear systems cannot be applied to analyze non-smooth oscillators. To understand the dynamic properties of non-smooth dynamical systems, the bilinear frequency approximation (BFA) method, a classical analytical tool, was developed to estimate the nonlinear natural frequency of bilinear oscillators by using the linear characteristics of underlying linear subsystems to construct a closed-form analytical solution. However, the traditional BFA method is only applicable to systems in which the gap size between the contacting element is zero, which limits the capability of this efficient analytical technique. In this paper, a generalized BFA method is proposed to extend the capability of the method to capture the nonlinear natural frequency of high-dimensional non-smooth systems where either gaps or prestress exist. The generalized BFA method employs the natural frequencies, mode shapes, and modal amplitudes of the linear subsystems as well as the gap size at the contacting degrees of freedom (DOF) to estimate the nonlinear natural frequencies at different amplitude levels. The backbone curve of non-smooth oscillators can be efficiently constructed using the proposed analytical method. The generalized BFA method is demonstrated on a single-DOF system, a three-DOF system, and a finite element model (FEM) of a cantilever plate with contacting elements. The results show that the proposed method provides an accurate approximation of the backbone curve for the first resonance, which has the most dominant vibration response, with very low computational effort, but a rather poor approximation for higher resonances. The proposed technique can be integrated with high-dimensional FEMs to analyze the nonlinear resonance of complex structures. The performance of the method is validated both numerically and experimentally in this work.

From period-doubling bursting to chaotic–periodic bursting in a modified Chua’s circuit

The aim of this paper is to reveal the dynamical mechanism of bursting oscillations in non-smooth dynamical systems, and the interest is focused on the effects of period-doubling bifurcation and chaotic attractor on the bursting oscillations. By modifying a fourth-order Chua’s circuit, a dynamical system possessing non-smooth switching manifold and multiple scale variables, is established. By taking the slow-varying variable as a bifurcation parameter, Subcritical non-smooth Hopf bifurcation, C-bifurcation, and period-doubling bifurcation are obtained in the fast subsystem. Further more, chaotic attractors, which are generated from a series of period-doubling bifurcations, are also observed. Based on the numerical simulations and bifurcation analysis of the fast subsystem, eight typical bursting patterns are obtained, and their dynamical mechanism is revealed. Study shows that period-2 limit cycle in the fast subsystem may cause the slow–fast system to oscillates in alternating small and large amplitudes; A complex spiking attractors made up of limit cycles and chaotic attractors may lead to spiking states changing from chaotic oscillations to periodic ones, and the trajectory may jump directly from a chaotic attractor to other stable one without undergoing any bifurcation of the fast subsystem.

Composite Poincaré mapping of double grazing in non-smooth dynamical systems: bifurcations and insights

Composite Poincaré mapping of double grazing in non-smooth dynamical systems: bifurcations and insights

Improving the monodromy matrix computation in pathfollowing schemes for nonsmooth dynamics

An enhanced pathfollowing strategy was recently proposed employing a modified Newton–Raphson solver accelerated by a low-cost routine via a Krylov subspace iterator. One of the key points of the numerical strategy is the computation of the Jacobian of the Poincaré map (i.e., the monodromy matrix), which, by serving as iteration matrix, governs both the accuracy and robustness of the whole strategy. In the context of nonsmooth nonlinear dynamic systems, this matrix can only be computed via numerical differentiation, which leads to subtractive cancellation errors.A new technique is here proposed by performing a central finite difference of the vector field associated with the ODE problems within the time-step integration to assemble the monodromy matrix. This technique allows us to bound the errors induced by the standard numerical differentiation.By means of a wide numerical campaign for meaningful multi-dof dynamical systems, the enhanced pathfollowing scheme is shown to yield a higher level of accuracy than the standard one, as well as higher robustness, with improved convergence properties for the whole numerical strategy.The C++ code implementing the proposed methodology is freely available at https://zenodo.org/record/7245478, and includes code for both the classic method as well as the new approaches for the computation of the Jacobian matrix.

A Nonsmooth Dynamical Systems Perspective on Accelerated Extensions of ADMM

Recently, there has been great interest in connections between continuous-time dynamical systems and optimization methods, notably in the context of accelerated methods for smooth and unconstrained problems. In this paper we extend this perspective to nonsmooth and constrained problems by obtaining differential inclusions associated to novel accelerated variants of the alternating direction method of multipliers (ADMM). Through a Lyapunov analysis, we derive rates of convergence for these dynamical systems in different settings that illustrate an interesting tradeoff between decaying versus constant damping strategies. We also obtain modified equations capturing fine-grained details of these methods, which have improved stability and preserve the leading order convergence rates. An extension to general nonlinear equality and inequality constraints in connection with singular perturbation theory is provided.

Numerical Analysis on Shimmying Wheels with Dry Friction Damper

The dynamics of the 1.5-degree-of-freedom model of towed wheel is investigated. Dry friction at the king pin is considered, leading to a non-smooth dynamical system. Beyond analytical and numerical linear stability analysis, the nonlinear vibrations are investigated by numerical bifurcation analysis with smoothing and by numerical simulations with event handling. The effect of dry friction at the king pin on the birth of separated periodic branches is presented on bifurcation diagrams. The presence of bistable parameter domains is also shown. The effect of smoothing is investigated by comparing bifurcation diagrams of the smoothed and the original non-smooth systems.

Existence and stability of the periodic orbits induced by grazing bifurcation in a cantilever beam system with single rigid impacting constraint

Grazing, which can induce many nonclassical bifurcations, is a special dynamic phenomenon in some non-smooth dynamical systems such as vibro-impact systems with clearance. In this paper, the existence and stability of the periodic orbits induced by the grazing bifurcation in a cantilever beam system with impacts are uncovered. Firstly, the Poincaré mapping of the system is obtained by adopting the discontinuous mapping method. Then, by means of shooting method, the periodic orbits are determined, followed by the acquisition of Jacobian matrix in the case of non-impact is obtained subsequently. In addition, for various impacting patterns, a combination of inhomogeneous equations and inequations is obtained to determine the existence of period orbits after grazing, with the stability criterion of the grazing-induced periodic orbits given. Numerical results verify the effectiveness of theoretical analysis. What is more, we also give a conjecture about the relationship between eigenvalues and the type of periodic orbits when eigenvalues are imaginary numbers.