{'fraction_non_alphanumeric': 0.04920220712734806, 'fraction_numerical': 0.01342602736997271, 'mean_word_length': 4.412614555256065, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Iterative Sequential Action Control for Stable, Model-Based Control of Nonlinear Systems

This paper presents iterative Sequential Action Control (iSAC), a receding horizon approach for control of nonlinear systems. The iSAC method has a closed-form open-loop solution, which is iteratively updated between time steps by introducing constant control values applied for short duration. Application of a contractive constraint on the cost is shown to lead to closed-loop asymptotic stability under mild assumptions. The effect of asymptotically decaying disturbances on system trajectories is also examined. To demonstrate the applicability of iSAC to a variety of systems and conditions, we employ five different systems, including a 13-dimensional quaternion-based quadrotor. Each system is tested in different scenarios, ranging from feasible and infeasible trajectory tracking, to setpoint stabilization, with or without the presence of external disturbances. Finally, limitations of this work are discussed.


I. INTRODUCTION
Model-predictive control (MPC) has been widely used over the past twenty years for control of linear and nonlinear systems [1], [2]. The rationale behind MPC is the following: at each time step, an L 2 or alternative variation of the cost function is locally optimized over time to obtain the openloop control as a function of time, only a small portion of which is actually applied to the system. The time horizon is then shifted, and the process is repeated at a later time step based on state feedback. Using this approach, it is clear that optimizing with respect to the entire open-loop control, although often effective, might not be the most efficient way to compute the control since most of the optimizer is typically discarded. In view of this observation, our previous work in [3]- [7] presented Sequential Action Control (SAC), a receding horizon approach for control of nonlinear systems, that exploits elements from hybrid systems theory [8], [9]. In contrast with the aforementioned MPC methods, the open-loop solution in SAC optimizes the needle variation [10] of the cost resulting in a single, constant magnitude action which does not optimize, but rather improves the cost function value across time steps. Since only a single action-obtained in closedform-is computed at each time step, control calculation is efficient. Our earlier work in [3]- [7] indicates that SAC can drive benchmark and challenging systems-including the cartpendulum, acrobot and pendubot in [3], hopping and humanoid locomotion in [5] and rotor vehicles with dynamics on Lie groups in [6]-close to a desired equilibrium. Nevertheless, it cannot achieve final stabilization and, as a result, switching to a locally stabilizing controller is often necessary [3], [6].
Authors are with the Neuroscience and Robotics Laboratory (N×R) at the Department of Mechanical Engineering, Northwestern University, Evanston, IL. Email: man7therakis@u.northwestern.edu, t-murphey@northwestern.edu. This paper presents iterative Sequential Action Control (iSAC), an extension of SAC that addresses one fundamental question that was left open in our previous work-how can SAC achieve consistent stabilizing behavior? We provide theoretical and practical results verifying that the iSAC method leads to closed-loop stability, while being efficient enough for application in situations with limited computational resources, e.g., mobile phones [5]. Simulation examples show that iSAC is consistently successful in a variety of control scenarios, ranging from benchmark inversion problems [11] to control of quadrotors performing complex three-dimensional maneuvers like flips. Unlike other system-specific methods (see Section VI), our approach leads to closed-loop stability (when the task is feasible) without modifications across examples; only the dynamic model that is being controlled is changed for each example.
In addition to the aforementioned points, the following novelties of iSAC distinguish this work from alternative MPC methods (see, e.g., [1], [2], [12]- [16] and references therein). In order to solve the open-loop problem, most MPC methods either employ nonlinear programming solvers (see [17] and [18] for a review) or solve a matrix Riccati differential equation as, for example, in [12]- [14]. On the contrary, the solution of the open-loop problem in iSAC has an analytic form which requires only forward simulations of the dynamics and a costate variable, which is computationally efficient and does not depend on nonlinear programming libraries. Moreover, control saturations can be incorporated without additional computational overhead. Finally, as opposed to many MPC alternatives that utilize discrete-time dynamics [1], [2], [14], [15], [19], iSAC uses continuous-time dynamics. As a result, variable-step integration may be used, which, combined with the previous points, significantly speeds-up the solution process.
To achieve closed-loop stability for iSAC, we apply a contractive constraint on the cost [20]- [23]. Contractive constraints have been widely used in the MPC literature to show closed-loop stability as an alternative to methods relying on a terminal (region) constraint [1], [2], [15], [16], [24], [25]. A disadvantage of these methods is that they require the computation of a terminal region, which has to be calculated separately for each system of interest. Many procedures have been proposed on how to find this region; see, e.g., [19], [24], [26]- [28] and references therein. On the contrary, the contractive constraint approach can be applied to a variety of examples without modification and is easily implementable in practice. In our case, a contractive constraint on the cost can be naturally integrated in iSAC, without additional work, by exploiting tools from hybrid systems theory, e.g., the mode insertion gradient (see [8], [9]), as will be explained in Section III. This paper is structured as follows: Section II provides a description of SAC; Section III describes iSAC and the features that differentiate this work from SAC. Sections IV and V analyze local and global stability properties of iSAC respectively. Section VI demonstrates applicability of iSAC to a variety of systems and control scenarios while conclusive remarks are given in Section VII. Finally, all proofs are presented in the Appendix.

II. PRELIMINARIES -SEQUENTIAL ACTION CONTROL
For convenience, we briefly summarize the work presented in [3]- [6]. SAC is a receding horizon method that enables real-time, closed-loop constrained control synthesis for broad and challenging classes of systems, including by following the cyclic diagram in Fig. 1a. We shall consider nonlinear systems with input constraints such thaṫ with u ∈ U and i.e., systems that can be nonlinear with respect to the state vector, x : R → X , but are assumed to be linear with respect to the control vector, u : R → U. The state will sometimes be denoted as t → x t; x(t i ), u(·) when we want to make explicit the dependence on the initial state (and time), and corresponding control signal. The SAC method uses objectives of the form with incremental cost l t, x(t) , terminal cost m t i + T, x(t i + T ) and time horizon T . Although (2) lacks a norm on control effort, this norm is included in the following steps as shown in (6). The SAC method calculates controls that optimally improve (not optimize) the objective (2) even if the objective is non-convex or unbounded. The following definition is necessary before introducing the open-loop problem solved in SAC.
Definition 1: An action A is defined by the triplet consisting of a constant control value, u A ∈ U, application duration, λ A ∈ R + and application time, τ A ∈ R, such that The open-loop problem P that is solved in each iteration of the receding horizon strategy can now be defined as follows Find action A such that J x(t; subject to  [17]. In order to solve (3), the SAC method follows four steps as illustrated above. (b) A sample, one-dimensional (m = 1) action schedule u * s (t) and the corresponding cost (8) used to calculate the application time τ A for the hypothetical time window [0, 4]s. Every point on u * s (t) corresponds to an action of infinitesimal duration, with value and application time as determined by the curve. Because u * s (t) is calculated from minimization of (6), all actions in u * s (t) would optimally improve (2) relative to only applying nominal control if applied for λ → 0 + . Each point on the mode insertion gradient curve approximates the change in cost (2) achievable by infinitesimal application of the corresponding action in u * s (t). By choosing the application time that minimizes this curve, we pick the action the generates the greatest cost reduction in the current time window.
From now on, subscript i will denote the i-th time step, starting from i = 0.
Definition 2: Nominal control u nom : R → U is continuous in t, and is used as a basis when calculating the open-loop solution. It is often u nom (·) ≡ 0 so that problem P outputs the optimal action relative to doing nothing (allowing the system to drift for a horizon into the future). Alternatively, u nom (·) may be an optimized feedforward or state-feedback controller, indicating potential shared control implementations. To evaluate SAC as a standalone controller, the nominal control throughout this paper is considered constant or zero. The system trajectory corresponding to application of nominal control will be denoted as x t; x(t i ), u nom (·) or x nom (·) for brevity.
The solution u * i (t) of problem P includes a switch to the calculated action (Fig. 2a) and thus, it is piecewise continuous in t. When applied to (1), u * i (t) generates a (discontinuous) switch of the same duration λ A in the dynamics from f t, x, u nom i (·) to f (t, x, u A ). The condition in (4) highlights a key feature of SAC, e.g., rather than optimize the objective (2), SAC actions improve (2) relative to only applying nominal control. As the receding horizon strategy progresses, P(t i , x i ) is solved from the current time t i using the measured state x i , and the calculated control u * i (t)-corresponding to x * i (t)-is applied for t s seconds (sampling time) to drive the system from x i to x t i + t s ; t i , x i , u * i (·) . The process is then repeated at the next sampling instance, i.e., t i ← t i + t s and i ← i + 1. The final result of the closed-loop receding horizon strategy is a sequence of actions, forming a piecewise continuous control signal u cl (t) with state response x cl (t).
Unlike alternative MPC methods [1], [2], [15], [16], [24], [25], the open-loop problem in SAC can be solved in closed form without employing nonlinear programming solvers (see [17] and [18] for a review). In order to solve (3), the SAC method follows four steps as illustrated in Fig. 1a. Beginning with prediction, the steps of the SAC process are described in the subsequent subsections.
Steps for solving the open-loop problem P The solution process assumes the following: Assumption 1: The control objective is to steer the state to the origin. This is not a restrictive assumption as most control scenarios (including trajectory tracking) can satisfy it with a change of coordinates.
Assumption 2: The dynamics f in (1) are continuous in u, piecewise continuous in t and continuously differentiable in x. Also, f is compact, and thus bounded, on compact sets X and U. Finally, the system is assumed to be controllable and f (·, 0, 0) = 0.
Assumption 3: The terminal cost m is positive definite, continuously differentiable in x and continuous in t. The incremental cost l(t, x) is continuous in t and continuously differentiable in x and l(·, 0) = 0. Also, there exists a continuous positive definite, radially unbounded function M : X → R + such that l(t, x) ≥ M (x) ∀t.
Assumption 4: The trajectory x * i (t) ∈ X corresponding to the solution u * i (t) ∈ U of P(t i , x i ) is absolutely continuous, and thus bounded, The open-loop problem P can then be solved by following the four steps presented below.
1) Predict: The SAC process begins by predicting the evolution of a system model from current state feedback. In this step, the system (1) is simulated from the current state x i and time t i , under u nom The sensitivity of (2) to the state variations along the predicted trajectory x nom i (·) is provided by an adjoint variable, ρ i : [t i , t i + T ] → R n , also simulated during the prediction step. The adjoint satisfieṡ where D i denotes derivative with respect to i-th argument. In short, the prediction phase completes upon simulation of the state using (1) and the adjoint system (5) under u nom i (·) control. The resulting trajectories x nom i (·), ρ i (·) will be used in the remaining three steps of the solution process.
2) Compute optimal action schedule u * s (·): In this step, we compute a schedule u * s : [t i , t i + T ] → R m which contains candidates for the next action, for all possible times that action may be taken. Specifically, u * s (·) contains candidate action values and their corresponding application times, assuming λ → 0 + for all (see Fig. 1b for a sample one-dimensional action schedule). The final u A and τ A will be selected from these candidates in step three of the solution process such that u A = u * s (τ A ), while a finite duration λ A will be selected in step four. The optimal action schedule u * s (·) is calculated by minimizing where the quantity dJ dλ (·), called mode insertion gradient (see [8], [9]), denotes the rate of change of the cost with respect to a switch of infinitesimal duration λ in the dynamics of the system. In this case, dJ dλ (·) shows how the cost will change if we introduce infinitesimal switches from f t, x nom Intuitively, minimization of (6) is driving dJ dλ (·) to a negative value α d ∈ R − . As a result, each switch/action value in u * s (t) is the single choice that optimally improves (2) (relative to only applying nominal control) if applied for λ → 0 + at its corresponding application time. The design parameter α d determines how much the cost is improved by each infinitesimal action in the schedule u * s (t). Based on the simulation of the dynamics (1), and (5) completed in the prediction step (Section II-1), minimization of (6) leads to the following closed-form expression for the optimal action schedule: where . The infinitesimal action schedule can then be directly saturated to satisfy any min/max control constraints of the form u min,k ≤ 0 ≤ u max,k ∀k ∈ {1, . . . , m} such that u * s ∈ U without additional computational overhead (see [3] for proof).
3) Determine application time τ A (and thus u A value): The SAC method optimizes a decision variable not normally included in control calculations-the choice of when to act. As opposed to always acting at the current time, i.e., τ A = t i , the application time of an action is allowed to take values in Recall that the curve u * s (·) provides the values and application times of possible infinitesimal actions that SAC could take at different times to optimally improve (2) from that time. In this step the SAC method chooses one of these actions to apply, i.e., chooses the application time τ A and thus an action Notice that the cost (8) is actually the mode insertion gradient evaluated at the optimal schedule u * s (·). Thus, minimization of (8) is equivalent to selecting the infinitesimal action from u * s (·) that will generate the greatest cost reduction relative to only applying nominal control. For a sample J t (τ ) see Fig. 1b. Finally, the selected application time must satisfy the following assumption: The role of this assumption is explained in the following step. 4) Determine control duration λ A : So far we have selected τ A and u A (see Definition 1) from a schedule of possible infinitesimal actions, u * s (·). The final step in synthesizing a SAC action is to choose how long to act, i.e., a finite control duration λ A , such that the condition in (4) is satisfied. The following assumption and proposition will facilitate the analysis in the sequel.
Proposition 1: For a choice of α d < 0 in (6), an infinitesimal control action u * s (τ ) that is selected according to Assumption 5 will result in dJ dλ (·) < 0 and thus, it can reduce (2). Proof: See Appendix. Proposition 1 proves that if α d < 0, then the infinitesimal u * s (τ A ) will lead to a negative dJ dλ (·). From [9], [29], there is a non-zero neighborhood around λ → 0 + where the mode insertion gradient models the change in cost in (4) to first order, and thus, a finite duration λ A exists that satisfies (4). In particular, for finite durations λ in this neighborhood we can write Thus, from Proposition 1 and (9), the condition in (4) is feasible. Then, a finite action duration λ A can be calculated by employing a line search process [29]. Starting with an initial duration λ 0 we simulate the effect of the control action using (1) and (2). If the simulated action satisfies (4), the duration is selected. If this is not the case, the duration is reduced and the process is repeated. By continuity, the final duration will produce a change in cost within tolerance of that predicted from (9). After computing the duration λ A , the control action A is fully specified (it has a value, an application time and a duration) and thus the solution u * i (t) of problem P has been determined. By iterating on this process (Section II-1 until Section II-4), SAC synthesizes piecewise continuous, constrained control laws for nonlinear systems. For more information about SAC, the reader is encouraged to consult [3], [5], [6].

III. FROM SAC TO ISAC
In this section we will describe the differences between SAC and iSAC, and the extensions introduced in iSAC that enable us to show closed-loop asymptotic stability. There are two main features introduced in iSAC; the first one is that the open-loop solution u * i (t) is iteratively updated as the time horizon progresses. The second is that the cost is decreasing from time step to time step. These two features are illustrated in Fig. 2. For convenience, the corresponding SAC behavior is also shown in the same figure.
A. Iterative update of u * i (t) Figure 2a illustrates how the iSAC method stores actions calculated at previous time steps and modifies u * i (t) by introducing a new action at every time step. On the contrary, SAC only keeps a single action in u * i (t), regardless of whether that action lies in the current application window [t i , t i + t s ] or not. In light of this, the following definition is necessary to distinguish between the nominal control (zero or constant in this paper) from Def. 2, and the default control that iSAC uses to calculate actions.
is piecewise continuous in t and is defined as with u def 0 (·) ≡ u nom 0 (·). In expression (10), As a result, actions calculated in iSAC improve the cost (2) with respect to applying default control u def i (·). In fact, as shown in Fig. 2, iSAC actions establish a decreasing behavior 1 for (2) from time step to time step. As will be explained in the following paragraph, this is accomplished by applying a contractive constraint on the cost, while the iterative nature of u * i (t) has a major role in the formulation of this constraint.

B. Contractive constraint on cost
In SAC, each action is calculated such that (4) is satisfied. However, in general, this does not ensure that the cost will follow a decreasing trend. The explanation lies in the fact that (4) and (9) only guarantee that the cost will improve relative to a nominal input in the current time step, but not necessarily across time steps. As an example, Fig. 2c shows a general scenario where condition (4) is met. In this case, the system will most likely not meet the desired specifications encoded in the cost, since the latter is not decreasing with time. Thus, even if (4) is satisfied, closed-loop stability cannot be guaranteed.
A solution to this problem is to modify condition (4) such that (2) sufficiently "contracts" between time steps. Specifically, instead of improving the cost relative to only applying nominal control, we could apply as a contractive constraint. Contractive constraints [20]- [23] have been widely used in the MPC literature to show closedloop stability as an alternative to methods relying on a terminal region constraint [1], [2], [15], [16], [24], [25]. Conditions similar to (11) (4). It can be observed that the behavior of the cost in SAC heavily depends on how J x(t; can take arbitrary values, there is no guarantee that the closed-loop system will be stable, even though condition (4) may be repeatedly satisfied.
in continuous or in discrete time, as an intermediate step used to prove closed-loop stability.
One problem that arises with the new condition (11) in SAC is that the quantities J x * i (·) and J x * i−1 (·) cannot be related through the mode insertion gradient, unlike, e.g., J x * i (·) and J x nom i (·) that appear in the original condition (4) and can be related through (9). Thus, Proposition 1 can no longer be used to ensure that (11) is feasible in SAC. On the contrary, in iSAC, we can use the iterative nature of u * i (t) to ensure feasibility of (11) through Proposition 1 for sufficiently small sampling time t s . In particular, in iSAC we can transform (11) to an equation that is similar to (4) and only includes the terms J x * i (·) and J x def i (·) , which can be related through the mode insertion gradient dJ dλ , as in (9). By definition, (10) Fig. 2b). We can then calculate Combining (12) and (11) and using (10) we get Thus, we were able to transform (11) into a sufficient decrease condition similar to (4), i.e., a condition that involves terms depending on the default and calculated open-loop trajectories. The following Proposition can be used in conjunction with Proposition 1 to ensure that the new condition (14) is (recursively) feasible.
Proposition 2: For α d < 0 there exists a sufficiently small sampling time t s such that the sufficient cost decrease required by the contractive constraint (14) is attainable.
Proof: See Appendix.
The contractive constraint (14) can be applied in the line search process that determines the duration of an action (Section II-4) in lieu of (4). By applying this condition at each time step, we can guarantee that the cost value will decrease across time steps, and using this fact we can prove closed-loop stability (Section V).

Open-loop problem in iSAC
Based on the above, the open-loop problem solved by iSAC at each time step is given by Find action A such that (14) is satisfied subject to Similar to the open-loop problem P in SAC, the solution u * i (t) of B can be obtained in closed form by following the four steps in Fig. 1a and Section II, without relying on nonlinear programming solvers. The only difference is that in the solution process of B, the superscripts nom are replaced by Proof: See Appendix.

IV. LOCAL STABILITY ANALYSIS
In this section we present local stability properties for iSAC, assuming quadratic cost (2). Note that this section is not a prerequisite for global stability analysis, which is presented in the following section. In other words, the results presented here are independent and not required to show global stability. The purpose of the section is to highlight local properties of iSAC that facilitate analysis of linear or linearized systems with quadratic cost.
Theorem 1: Assume that the cost (2) is quadratic with symmetric, positive semidefinite Q and P 1 . Then, there exist neighborhoods around the equilibrium and final time, N (0) and N (t i + T ) respectively, where (7) simplifies to a linear feedback regulator where P (t) is a n × n matrix that can be calculated from a Riccati equation.
Proof: The proof can be found in the Appendix. Similarly, one can also show Theorem 1 yields a feedback expression in error coordinates for trajectory tracking problems. The following Corollary shows that the quantity P (t) is calculated by solving the same Riccati differential equation that runs backward in time as in the LQR solution [30].
Corollary 1: The optimal action schedule (17) is locally equivalent to the LQR solution u * Thus, (17) also minimizes (18). Proof: The proof can be found in the Appendix. Remark 1: Notice, that since the schedule (17) minimizes (18), if (17) is directly applied to the system as a sequence of optimal infinitesimal actions, i.e., set u * i (t) = u * s (t), the closed loop system will inherit the properties of an LQR response. In this case, one can use LQR theory [30] to prove existence of a Lyapunov function that guarantees local closedloop stability for iSAC. Additionally, there is no need to search for an application time and duration in this case as (17) models a sequence of optimal infinitesimal actions. Instead, one has to calculate P (t) by solving a Riccati differential equation. Naturally, the contractive constraint (14) does not apply here either since (17) minimizes (18) as opposed to "contracting" it. Switching to this linear, duration-independent equivalent to iSAC near the equilibrium can be useful, e.g., in cases where the duration λ required to satisfy (14) is shorter than the minimum tolerance imposed by numerical software or hardware. Sums-of-Squares (SOS) and the S-procedure [28], [31] can pre-compute regions of attraction for (17) and so can determine when the switch to this duration-independent version of iSAC should occur.

V. GLOBAL STABILITY ANALYSIS
In this section, we provide global stability results for iSAC and discuss how iSAC can be used under a special case of disturbances and to track infeasible trajectories, i.e., trajectories that do not satisfy the dynamics (1).

A. Nominal Case
In this section we are considering nominal stability (not robust stability), i.e., the trajectories of the plant are assumed to coincide with the trajectories predicted by the model. Recall that the iSAC method does not optimize the cost (2) but rather improves it across time steps. It is often the case that minimization of the objective function in the open-loop optimization problem is not required to achieve stability of a model predictive controller. Sometimes, a decrease in the cost at every iteration is sufficient to guarantee stability [19], [24], [32], [33]. In view of these observations, we now present the following Theorem: Theorem 2 (Asymptotic Stability): For sufficiently small sampling time t s , the closed-loop system resulting from the receding horizon strategy is asymptotically stable in the sense that ||x cl (t)|| → 0 as t → ∞.
Proof: See Appendix. Theorem 2 does not imply Lyapunov stability, but rather establishes the usual notion of attractiveness. The importance of t s in establishing closed-loop stability is also underlined; if the assumptions of the Theorem hold, by appropriate selection of t s , we can ensure that the contractive constraint is feasible, and thus, asymptotic stability.

Cost as a Lyapunov function
Now that asymptotic stability has been proved, we will show that, under certain assumptions, the objective (2) is a Lyapunov function for the closed-loop system that decreases at intervals of prediction horizons.
Theorem 3 (Objective Function as a Lyapunov Function): Assume that f is bounded such that ||f || ≤ ξ||x|| for some finite constant ξ > 0, and that the cost (2) is of the form in (16), rewritten as Then, for sufficiently small t s , (19) is a discrete Lyapunov function for the closed-loop system. Proof: The proof can be found in the Appendix. Corollary 2: For the assumptions in Theorem 3, if t s → 0, then (19) is a continuous Lyapunov function for the closedloop system, implying asymptotic stability in the classical sense, i.e., Lyapunov stability and attractiveness.
Proof: The proof can be found in the Appendix.

B. Asymptotically Decaying Disturbances
This section provides stability results for a special case of disturbances. Here, we assume that the plant is described by the following differential equations (instead of (1)) x p = f (t, x p , u) + η(t) = g(t, x p ) + h(t, x p ) u + η(t), (20) where x p (t) represents the trajectory resulting from applying control u and the effect of an unknown additive disturbance η(t). The disturbance is bounded and asymptotically decaying; in particular, we assume the following: Furthermore, η(t) is asymptotically decaying, i.e., δ i → 0 as i → ∞. Closed-loop stability under η(t) is provided by the following theorem.
Theorem 4 (Asymptotic Stability under Asymptotically Decaying Disturbances): For sufficiently small t s , the closedloop system under the disturbance η(t) is asymptotically stable in the sense that ||x cl (t)|| → 0 as t → ∞.
Proof: The proof can be found in the Appendix. As seen in the proof, for as long as the disturbance is acting, the cost is allowed to increase, but the amount of increase is kept sufficiently bounded by iSAC. As time progresses, the disturbance attenuates and the contractive constraint becomes feasible again.

C. Remarks on tracking of non-admissible trajectories
When designing a control task, there are often cases where asymptotic stability cannot be achieved. For example, in trajectory tracking, it is not always easy to identify trajectories that satisfy the dynamics of the system of interest, unless there is special structure one can exploit, e.g., differential flatness [34]. The control goal in this case is to ensure that the deviation of the generated (feasible) open-loop trajectory from the desired one is bounded, i.e., ||x * Unlike the decaying disturbance case, the upper bound ∆ i does not have any structure that can be exploited. Thus there is no guarantee that x cl (t) → x d (t) as t → ∞.
Similar to the disturbance case, these observations stem from the fact that the feasibility of the contractive constraint, and thus of the open-loop problem, can no longer be ensured when tracking an infeasible trajectory. In these cases, in order to keep the error from the desired behavior bounded, one can set α d to more negative values such that the control generation is more "aggressive". By doing so, the selected action leads to smaller cost increase (at a higher control cost) and iSAC generates a feasible response that is close to the desired infeasible objective, as shown in the simulations that follow.

VI. SIMULATION RESULTS
In this section we demonstrate the flexibility and versatility of iSAC by applying the theoretical results presented in the previous sections to a variety of challenging systems. It is important to note that unlike other system-specific methods mentioned in this section, our approach requires no modifications across examples; only the dynamic model that is being controlled and the relevant parameters of iSAC are changed for each example. Also, all the examples that will be presented here run much faster than real time, indicating applicability to hardware implementations. For better visualization of the results, the reader is encouraged to visit https://vimeo.com/219702474. The control objective is encoded using a quadratic cost as in (16)

throughout this section (Assumption 3 is thus satisfied).
To demonstrate Theorem 1, the control generation in the first example switches to the state feedback version of iSAC in (17) near the equilibrium. The remaining examples use iSAC without switching to state feedback and also include cases of infeasible tracking tasks and disturbance rejection. The sampling time t s for each example was appropriately selected to achieve closed-loop stability (Theorems 2 and 4) when feasible. Finally, a summary of iSAC-specific parameters for each example is given in Table I, while Fig. 3 shows the resulting trajectories for the examples described below. Our method is shown to be successful in these scenarios, leading to asymptotic stabilization when feasible.

A. Acrobot
The dynamics of the system and the parameters used in the simulations that follow can be found in [35]. The state vector for the system includes the two joint angles and their Swing-up problem: The swing-up problem for the acrobot is a benchmark problem in control theory, where the system starts from the downward equilibrium and the control objective is stabilization at the upward equilibrium [36], [37]. The desired setpoint was thus x d = [0, 0, 0, 0]. To illustrate the durationindependent version of iSAC presented in Section IV, control calculation was based on Theorem 1 after system trajectories were within a pre-defined region of the upward equilibrium, found also in [35]. The cost weights were selected as Q = Diag[1000 , 0 , 500 , 11], and R = 0.01. Figures 3 and 4 show the resulting trajectories, closed-loop control and cost generated by iSAC. In contrast with other methods that switch between separately derived controllers [38], iSAC pumps energy into the system by appropriate excitation, without explicitly encoding this behavior in the objective. Finally, many methods (see [36], [39]) rely on energy-based objectives to solve the swing-up problem for the acrobot and pendulum systems, primarily because statetracking objectives result in many local minima that prevent convergence to desirable trajectories. It is noteworthy that, like SAC, iSAC is able to invert the system using a statetracking objective while bypassing local minima. For a detailed comparison between SAC, sequential quadratic programming (SQP) [40], and iLQG [14] the reader is encouraged to consult [3] and [6].

B. Cart-pendulum
The dynamics of the system and the parameters used in the simulations that follow can be found in [35]. The state vector for the system includes the angle and angular velocity of the pendulum and the horizontal position and velocity of the cart, [θ,θ, x c ,ẋ c ]. The control input is the horizontal acceleration of the cart.
2) Feasible trajectory tracking: The desired trajectory in this example was selected as θ d (t) = π 5 sin t and the cost weights were Q = Diag[5000 , 0 , 0 , 0], P = Diag[500 , 0 , 0 , 0], R = 0.3. The desired trajectory was feasible given that no constraints were applied to the cart position and velocity. As the plot in Fig. 3 illustrates, tracking was very accurate, making it hard to distinguish between the desired and actual trajectory followed by the system. The corresponding cost is shown in Fig. 4.
As mentioned in the acrobot example, it is worth noting that iSAC is able to solve the inversion problem by pumping energy into the system without relying on separate controllers. In addition, this was accomplished using the state tracking objective, as opposed to the energy objective commonly employed (see for example [11]), bypassing local minima that cause solutions to converge to undesired trajectories in state tracking costs.

C. Car-like vehicle
The dynamics used in the simulations arė where x c , y c denote the position of the car, θ is the angle with respect to the horizontal axis and v is the forward velocity. The iSAC method directly controls the acceleration and the angular velocity of the car. This nonholonomic system violates Brockett's necessary condition for smooth or even continuous stabilization [42], which makes the control design problem challenging. Since iSAC automatically generates a discontinuous control law, it is not subject to this condition.
2) Feasible trajectory tracking under disturbance: The desired trajectory in this example was chosen as x d c (t), y d c (t), θ d (t) = (t, sin t, tan −1 cos t) and the cost weights were Q = Diag[100 , 100 , 0 , 10] and R = Diag[0.1 , 0.1]. An external disturbance acted instantaneously on the system on four occasions, i.e., at t = 0, 3, 6, and 9s, each perturbing v and θ by +3 m/s and + π 4 rad respectively. The effect of the disturbance can be seen in the corresponding plot in Fig. 3; while the car was tracking the desired trajectory, the disturbance pushed the system away every time it acted on the system. After the last disturbance attenuated the control successfully steered the system back to the desired trajectory. The corresponding cost is in Fig. 4.
A great deal of work on nonholonomically constrained carlike models was completed in the 1990s using nonsmooth or time-varying control laws aiming to overcome stabilizability limitations traditional techniques [43]- [45]. Alternative methods that achieve similar results in the parking problem include fuzzy controllers as in [46] and dynamic feedback linearization [47]. Other approaches that intrinsically lead to piecewise continuous controls like iSAC include MPC as in [20], [48]. The former paper utilizes contractive constraints on the state and, although successful, results in a straight line path between the initial and final configuration. The latter uses the aforementioned terminal region constraints to achieve stabilization, but, unlike iSAC, the method is designed specifically for the car-vehicle system.

D. Planar vertical take-off landing (PVTOL) aircraft
The dynamic model used in this example is given in [49]. The state vector for the system includes the position of the aircraft, the angle with respect to the horizontal axis and the corresponding velocities, i.e., [x a ,ẋ a , y a ,ẏ a , θ,θ]. The control inputs of the system are thrust (directed out the bottom of the aircraft) and the rolling moment. Lastly, the parameter that appears in the model was chosen as = 0.3 and = 0.05 for the two examples that follow respectively. Note that couples the rolling moment input with the lateral acceleration of the aircraft; a positive value > 0 means that applying a (positive) moment to roll left produces an acceleration to the right, making this a non-minimum phase system.

1) Gust control:
The desired setpoint in this example was x d = [1, 0, 0.5, 0, 0, 0] and the cost weights were selected as Q = Diag[15 , 3 , 15 , 3 , 3 , 0] and R = Diag[0.1 , 0.1]. An external disturbance acted instantaneously on the system on five occasions, i.e., at t = 0, 2, 4, 6, and 8s, each perturbingẋ a by +0.1 m/s. The behavior of the system was similar to the one observed in the car-like vehicle; while the control was steering the system towards the desired setpoint, the disturbance drove it away. After the last disturbance decayed and the setpoint was finally reached.
2) Infeasible figure eight tracking: The desired trajectory for this example was a planar figure eight. To track this trajectory, the aircraft needs make aggressive maneuvers; it is therefore impossible to keep the angle zero and track the desired position at the same time which makes the selected trajectory infeasible. Nevertheless, as seen in Fig. 3 Figures 3 and 4 show the resulting trajectories and cost. The control of PVTOL aircrafts has been extensively studied using, e.g., gain-scheduling [50], robust control [51] and input-output linearization [49]. In the latter case, the system is decoupled into simpler subsystems to facilitate the control process (a similar approach is commonly followed in quadrotor control). It is important to note that iSAC is able to successfully control the system using the full dynamics and without relying on separate controllers, unlike, e.g. in [50]. Finally, other techniques that have been applied on VTOL aircrafts and unmanned aerial vehicles (UAVs) in general include sliding mode control and backstepping as explained in the following example [52].

E. Quadrotor
For this example, a quaternion-based model is used because it leads to polynomial, singularity-free dynamical equations for the quadrotor. The model can be found in [53] and the model parameters we used are the same as in [6]. The 13-dimensional state vector consists of the position and velocity of the center of mass, the angular orientation with respect to the inertial frame expressed in quaternions and the angular velocity expressed in the body frame, i.e., [x q ,ẋ q , y q ,ẏ q , z q ,ż q , q 0 , q 1 , q 2 , q 3 , p, q, r]. The control inputs are the squared angular velocities of the rotors, which are converted to upward thrust for each rotor when multiplied by an appropriately dimensioned constant. 2) Infeasible three dimensional figure eight tracking: In this example, the quadrotor starts in an upsidedown position and is set to track a three-dimensional figure eight while keeping the quaternion states at In other words, the quadrotor is requested to track an aggressive 3-dimensional maneuver without changing its angular orientation, which is infeasible. The cost weights were selected as Q = Diag[10 , 10 , 10 , 0 , 0 , 0 , 1000 , 1000 , 1000 , 1000 , 0 , 0 , 0] and R = Diag[0.1 , 0.1 , 0.1 , 0.1]. As seen in Fig. 3, at the beginning of the simulation the quadrotor reverts its angular orientation and then it starts tracking the desired trajectory keeping the error small. This behavior is automatically generated by our method. The corresponding cost is in Fig. 4.
Many different approaches have been followed for control of quadrotors over the years. Linear control methods like LQ and PID synthesis [54], [55] are popular for their simplicity. They utilize decoupled or simplified dynamics and also exploit differential flatness [34] in order to track feasible trajectories. For example, in [56] a quadrotor flip is implemented using a dedicated attitude controller. In contrast iSAC controls the full nonlinear dynamics of the system, and has good performance even with infeasible trajectories as shown in the simulation examples. No additional steps are required to apply our method to this system, unlike, e.g., feedback linearization [57] or sliding mode control [52]. The simplicity and versatility of iSAC is one of the main points that this section demonstrates; the same exact approach is being applied to a variety of challenging systems in real time, by replacing only the dynamic model that is being controlled.

F. Discussion and Limitations
In general, iSAC performs well in the variety of scenarios presented. When the task is feasible, trajectories are asymptotically stable and the cost is generally decreasing from time step to time step as a result of applying condition (11). Nevertheless, near the desired equilibrium, the behavior of the cost sometimes deviates from what one would expect, even when the control task is feasible. As an example, in the rightmost panel of Fig. 4 the cost is not continuously decreasing near the equilibrium, although, in both of these examples, the system is effectively stabilized when the unexpected behavior occurs (the cost value is less than 0.01). There are three ways to explain this; first, in order to speed up the solution process, the backtracking process (final step in the solution of the open-loop problem) in all examples presented in this paper was stopped after 10 iterations. While for the majority of cases 10 iterations were sufficient to satisfy the sufficient descent condition (11), on some occasions more iterations were necessary. As a result, a small increase was observed in the cost in these time steps. Another factor to consider is numerical tolerance. For example, from the simplified model (9), depending on the value of the mode insertion gradient, it could be the case that the required λ value that achieves the sufficient decrease descent is smaller than the minimum tolerance of numerical integrators. This issue appears often near the equilibrium where systems are more sensitive to the duration of an action. One way to deal with this issue is to switch to the duration-independent, state feedback version of iSAC near the equilibrium (Theorem 1). Finally, from Proposition 2, it is possible that the selected sampling time t s is not sufficiently small to achieve (recursive) feasibility of (14). In this case, selecting a smaller t s would likely solve the problem.

VII. CONCLUSION
In this paper we presented iSAC, a model-predictive approach for control of nonlinear systems. As the time horizon progresses, our method sequences together optimal actions and synthesizes piecewise continuous control laws. Some key characteristics of iSAC include: a) there is no need to rely on nonlinear programming solvers since there is an analytic solution to the open-loop problem, b) it is iterative in the sense that at each time step it updates the open-loop solution of the previous time steps by adding one action and c) it utilizes continuous dynamics while incorporating control constraints without additional overhead. Due to these points, iSAC leads to computationally efficient solutions. To establish closed-loop stability, we applied a contractive constraint on the cost. Compared to methods relying on terminal region constraints, the contractive constraint alleviates the need to calculate a terminal region. We also investigated different control scenarios ranging  I   ISAC PARAMETER VALUES USED IN SIMULATIONS   Example  VI-A  VI-B1  VI-B2  VI-C1  VI-C2  VI-D1  VI-D2  VI-E1 VI- As explained in Section VI-F, numerical issues may arise when the system trajectories are near the desired equilibrium (see e.g., gray area in the plots). As a result, the contractive constraint may not be satisfied momentarily leading to a small cost increase. While this may not be a problem since the system is effectively already stabilized at that point (compare corresponding trajectories from Fig. 3), one way to deal with this issue is to switch to the duration-independent, state feedback version of iSAC near the equilibrium (Theorem 1) as demonstrated by the acrobot swing-up panel. Finally, columns three and four show the cost behavior when acting under the presence of an external bounded disturbance and tracking infeasible trajectories respectively. For the former, the unmodeled disturbance leads to a temporary cost increase, but our method responds accordingly as the disturbance attenuates. For the latter, as explained in Section V-C, the contractive constraint is no longer feasible which explains why the cost is not consistently decreasing. Nevertheless, iSAC keeps the deviation from the desired trajectory, e.g., the cost, bounded.
from feasible and infeasible trajectory tracking to set point stabilization with or without external disturbances. Finally, we presented simulation examples using five different systems to demonstrate the applicability and flexibility of our method.

B. Proof of Proposition 2
From Proposition 1 and (9), we know that J x * i (·) − J x def i (·) = ∆J < 0 for each λ in a neighborhood around λ → 0 + if α d < 0. For the special case of t s → 0, it follows directly from Proposition 1 that (14) is always feasible since C ≡ 0. For the general case of nonzero sampling time, by continuity assumptions for m and l, C in (14) is continuous with respect to t s . Thus, it follows that there exists a sufficiently small t s such that ∆J < C.

C. Proof of Proposition 3
The open-loop problem B is solved by following the same four sequential steps as in the open-loop problem P in Section II. The only difference is that in the solution process of B, the superscripts nom are replaced by def. To show that u * i (t) exists, we will show that each of the four solution steps has a solution: 1) The first step in Section II-1 involves calculating the solutions to (1) and (5). Since the default control is in general discontinuous, if solutions are interpreted as sample and hold (CLSS) solutions (see [58]), existence follows directly from Assumptions 2, 3.
2) In the second step our method in II-2 calculates the optimal action schedule u * s by minimizing (6). Because (6) is convex with a continuous first variation from Assumptions 1-3, solutions (7) exist and are unique, which is also both necessary and sufficient for global optimality of (6).
3) In the third step in II-3 the process selects τ A and u A by minimizing (8). Because τ A ∈ [t i , t i + T ] and (8) is in general piecewise continuous (and thus bounded), a solution to this one-dimensional problem exists. 4) Finally, from Proposition 2, the backtracking process in Section II-4 is guaranteed to find a duration that satisfies condition (14) for sufficiently small t s .
Since all four subproblems have solutions, the open-loop solution u * i (t) exists.

D. Proof of Theorem 1
At final time, ρ(t i + T ) = P 1 x(t i + T ) from (5) and (16). Due to continuity, the following linear relationship must exist for a nonzero neighborhood of the final time, N (t i + T ), Ignoring the subscripts and time dependence for brevity, the schedule (7) can be formulated using (23) as For x ∈ N (0) the quadratic terms in the relationship above go to zero faster than the linear terms, and thus u * s (t) converges to the linear feedback regulator u * s (t) = u def + R −1 h(t, 0) T P (t)x(t)α d .

E. Proof of Corollary 1
This can be shown by Direct application of the Maximum Principle to (18).

F. Proof of Theorem 2
The proof has two parts; Lemma 1 shows that the integral t t0 M x(s) ds is bounded for t → ∞ (see Assumption 3). The latter is used in conjunction with a well-known lemma found, e.g., in [59], [60], to prove asymptotic convergence in the second part of the proof.
Consider the horizon interval [t i , t i + T ]. Let u * i (t) be the solution to the open-loop problem B(t i , x i ) and x * i (t) the corresponding state trajectory. Clearly, V t i , t i + T, x * i (·) = J x * i (·) . Then, for t ∈ [t i , t i + T ] V t, t i + T, x * i (·) Based on the above, we now present the following lemma. Lemma 1: For small t s , all t ∈ [t i , t i + T ] and all i ∈ N t t0 M x(s) ds ≤ V t 0 , t 0 + T, x * 0 (·) with Proof: From (11) and Assumption 3 we get that − V t i , t i + T, x * i (·) . (28) Also, from (25) we have − V t, t i + T, x * i (·) . (29) Adding (28) and (29) leads to (26) and Lemma 1 is proved.
From Lemma 1, because V t, t i + T, x * i (·) ≥ 0 and M is positive definite, we can deduce that t t0 M x(s) ds is bounded for t → ∞. We also have that x * i (·), and thus x cl (·), are bounded and from the properties of f ,ẋ * i (·) andẋ cl (·) are bounded as well. These facts combine with the following well-known lemma to prove asymptotic convergence.