# Approximate Prediction-Based Control Method for Nonlinear Oscillatory Systems with Applications to Chaotic Systems

Authors: Thiago P. Chagas; Pierre-Alexandre Bliman; Karl H. Kienitz
Journal: Journal of Control Science and Engineering
 (2018)

---

## Abstract

The approximate Prediction-Based Control method (aPBC) is the continuous-time version of the well-known Prediction-Based Chaos Control method applied to stabilize periodic orbits of nonlinear dynamical systems. The method is based on estimating future states of the free system response of continuous-time systems using the solution from the Runge-Kutta implicit method in real time. Some aspects of aPBC are evaluated in the present work, particularly its robustness to low future states estimation precision is exemplified.


---

## Body

## 1. Introduction

Oscillatory systems are typical in many problems of engineering, biological sciences, physics, economy, and other areas [1–7]. In general, oscillations shall be damped in order to reduce amplitude or stopped avoiding damage, reducing costs, increasing precision, and others [1–4]. These oscillations are mainly divided into periodic and aperiodic. Chaotic aperiodic oscillations are related to unpredictability, disorder, and instability while periodic oscillations are associated with order [8, 9].Chaotic sets are composed of an infinite number of unstable periodic orbits (UPOs) [10, 11] and the stabilization of one of these orbits leads to periodic oscillation, possibly, reducing amplitude. Different chaos control methods have been developed aiming at stabilizing these orbits with low control effort using the main characteristics of chaos [9].One classical chaos control method is the Delayed Feedback Control (DFC) proposed by Pyragas (1992) [12]. This is a state-feedback method whose control signal is computed through the difference between delayed and current time measured system states. DFC was initially proposed for continuous-time and has many applications on discrete-time systems [13]. Nevertheless, DFC has a well-known limitation proofed for discrete-time systems, the odd-number limitation [14–17], which means that DFC cannot stabilize orbits with an odd number of Floquet multipliers larger than +1. This limitation is questioned in literature for continuous-time systems with counter example [18–20].Different modifications on DFC have been proposed to overcome the odd-number limitation [13] and one of the most interesting modifications is the Prediction-Based Control (PBC) [21] developed for discrete-time systems. Instead of the delayed state, it uses the future state one period of the target UPO ahead computed along the trajectories of the free system response as reference for the control signal. Developments of PBC are proposed always for discrete-time systems; for example, in [22] a method is proposed for tuning control gain without previous knowledge about the target period-p UPO that leads to fast convergence rate of trajectories to the stabilized periodic orbit; in [23] it is proved that period-2 UPOs can be stabilized by a pulsed control signal reducing the required control effort for stabilization; and in [24] equilibrium points (or period-1 UPOs) are stabilized in the presence of multiplicative or additive noise.The literature provides some results on a PBC-like strategy applied to continuous-time systems. In [25] it is proposed to use the system Jacobian matrix at each point of the trajectory, instead of system future states, for stabilization of equilibrium points. The results showed that this strategy does not lead, necessarily, to stabilization of free system equilibrium points or UPOs and control effort does not vanish in steady state, two main advantages of chaos control proposed by Ott et al. [26]. In [27], the strategy proposed in [25] is incremented using Neural Networks also for stabilization of equilibrium points.In fact, PBC has a practical limitation for continuous-time systems if it is applied as proposed for discrete-time systems [21], the need of future state values in real time. In [28], we proposed the approximate Prediction-Based Control (aPBC) with a methodology based on the implicit Runge-Kutta method and state estimation applied to predict future states for the free system in real time based on system model. The authors claim that aPBC is the continuous-time applicable version of PBC because it uses prediction and stabilizes free system UPOs, ideally, vanishing steady state control effort. On the contrary, predicting future states has some drawbacks because the proposed future state prediction scheme requires an increase in the closed-loop system order and consequently computational power for numerical integration. Moreover, it is also subjected to prediction model and real system mismatch and integration method precision for application.Continuing the aPBC development and study, one of the main questions that remained in [28] is about the robustness of the method for low future state estimation precision. The forced van der Pol (vdP) oscillator (nonautonomous) and Rössler system (autonomous) are used in the present work for numerical examples on the aPBC robustness and trade-off between estimation accuracy and computational cost. We present a methodology and performance indexes that aid finding this trade-off and show results evidencing that it is possible to find a lower bound for precision. Besides that, this work presents an UPO of the forced van der Pol oscillator with one of the Floquet multipliers larger than +1 that is stabilized by the aPBC and could not be stabilized by the DFC with constant control gain following the procedures proposed in [18–20]. The optimization procedure proposed in [29] was also applied to the DFC without success. This is evidence that the aPBC maintain the advantage of discrete-time PBC of not being sensitive to the odd-number limitation.The paper is divided as follows. The aPBC is reviewed in Section3, generalizing its proposition. In Section 4, aPBC application is presented using orthogonal collocation method [30, 31] as implicit Runge-Kutta method. In Section 5, the forced van der Pol (vdP) oscillator and Rössler system and their chaotic behavior are presented. Both systems are used for numerical examples on the aPBC robustness and trade-off between estimation accuracy and computational cost in Section 6. The example of an UPO of the forced vdP oscillator that could not be stabilized by DFC and is stabilized by aPBC is shown in Section 6.1.4.


## 2. Problem Statement

Consider the following continuous-time dynamical system:(1)x˙t=ft,xt,ut,x0given,where t∈R+,  x:R+→Rn,  u:R+→Rm,  n,m∈N, and f:R+×Rn×Rm→Rn is a T-periodic function with respect to time t; that is, by definition(2)ft+T,x,u=ft,x,u,∀t∈R+,∀x∈Rn,∀u∈Rm.Moreover, we assume the existence of aT-periodic solution x∗(t) to the free system (1), which is the system obtained by setting u(t)=0,  t≥0. In other words,(3)x∗t+T=x∗t,∀t∈R+,x˙∗t=ft,x∗t,0m,∀t∈R.We assume that this periodic solution is unstable and can be stabilized by the continuous-time version of the PBC with feedback law defined as(4)ut,xt=Kt,xtφt+T,t,xt,0-xt,where u:R+×Rn→Rm is the control signal, K:R+×Rn→Rm×n is the control gain, φ(t1,t0,x,0) is the value at time t1 of the state of (1) with x(t0)=x and u(t)=0,  t0≤t≤t1. In other words, φ(t1,t0,x,0) is the value at time t1 of the state along the trajectory departing from x at time t0 of the free system (u(t)≡0).For an appropriate control gainK(t,x(t)),  x∗ is a stable solution of the closed-loop system(5)x˙t=ft,xt,ut,xtand control signal u ideally verifies(6)ut,x∗t=0t≥0.The general formulation of the continuous-time PBC method can be written as partial differential equation where the solution of system (5) with control signal (4) is the solution of the PDE (7a) and (7b). Consider(7a)∂Xt,0∂t=ft,Xt,0,Kt,Xt,0Xt,T-Xt,0,t≥0(7b)∂Xt,Ts∂s=Tft+Ts,Xt,Ts,0,t≥0,s∈0,1(8)X0,0=x0. The function X:R+×0,T→Rn is such that X(t,0)=x(t) and X(t,T)=φ(t+T,t,x(t),0).Condition (6) ensures zero control effort when the trajectory is on the unstable periodic solution x∗ of the free system. However, applying u(t,x(t)) in real systems is based on the exact prediction and real-time computation of future states φ(t+T,t,x(t),0). Due to these practical constrains, alternatively, we are also interested in stabilizing an orbitclose to the UPO x∗ of (1) resulting in a low control effort. This task is solved using the aPBC reviewed in the sequel.


## 3. Principles of the Approximate Prediction-Based Control Method (aPBC)

Computing the future state of the prediction termφ(t+T,t,x(t),0) requires solving the free system ODE at each time t, from time t to t+T. This cannot be done exactly in real time and its solution is the basis of the aPBC. This task is divided into two steps, the first one consists in approximating the solution of φ(t+T,t,x(t),0) by an implicit Runge-Kutta (R-K) ODE integration method, and the second is expressing the solution of the R-K method as a state observer that can be integrated by any explicit method in real time.


### 3.1. Approximation of the Prediction Term: 1st Step

The first step consists in approximating the solution of (7b) by an implicit R-K ODE integration method [32], in order to estimate the prediction term, that is, the terminal value,(9)Xt,T=x+T∫01ft+Ts,Xt,Ts,0ds,xgiven.In order to estimateX(t,T) given by (9), the state transition map of the free system φ(t2,t1,x,0) is first approximated by the operator z defined by(10a)zt+T,t,x=x+T∑i=1Ncilit(10b)lit=ft+Tsi,x+T∑j=1Naijljt,0,where lj:R→Rn and aij,ci∈R are weights chosen according to the implicit method used [32]. The approximation z(t+T,t,x(t)) of X(t,Ts),  s∈[0,1], is calculated at the discretization points s=si,  i=1,…,N. N shall be chosen accordingly such that the sums in (10a) and (10b) lead to an adequate approximation of the integral in (9) for the given f(t,x,u).For simplicity, (10b) is written in the vector form (11)Lt=FTt,x,Lt,where (12)Lt=l1t⋮lNt∈RnN,∀t≥0, and FT:R×Rn×RnN→RnN is defined by (13)FTt,x,L=ft+Ts1,x+∑j=1Na1jlj,0⋮ft+TsN,x+∑j=1NaNjlj,0,∀x∈Rn,∀L∈RnN.In order to computez(t+T,t,x) through (10a), it is necessary to solve the algebraic system of (11) with unknown L(t)∈RN.Writing(14)C=c1⋯cN and closing (5) by (15)ut,xt=Kt,xtzt+T,t,xt-xt=TKt,xtCLt yield the differential algebraic equation (DAE)(16a)x˙t=ft,xt,TKt,xtCLt,x0given(16b)Lt=FTt,xt,Lt.The real-time solution of the DAE (16a) and (16b) requires the computation of L(t), its algebraic term, at each time t. We therefore introduce an observer equation in the sequel, to transform the controlled system into a system of ODEs.


### 3.2. Approximation of the Prediction Term: 2nd Step

We now approximate (16b) by solving the nN-dimensional ODE (17) whose solution L^(t) is an estimation of L(t). The initial value L^(0) is intended to be (precisely) computed off-line to provide good tracking quality for L(t).(17)ddtL^t-FTt,xt,L^t+koL^t-FTt,xt,L^t=0,L^0=L0.The scalar gainko is chosen positive in order that L^(t) of (17) tends asymptotically towards the solution L(t) of (16b) when t→+∞ and typically in such a way that the estimator dynamics are faster than the controlled system dynamics. If the evolution of L^(t) may indeed be chosen in order to fulfill (17), convergence does occur.From (17), we obtain(18)L^˙t-∂1FTt,xt,L^t+∂2FTt,xt,L^t1N⊗ft,xt,0+∂3FTt,xt,L^tTA⊗InL^˙t+koL^t-FTt,xt,L^t=0,where A=(aij) and ∂i is the partial derivative with respect to the ith variable. In is the n×n identity matrix and 1N is the column vector of dimension N with all elements equal to 1. ⊗ is the Kronecker product.By determiningL^˙(t) from (18) we get (19)L^˙t=InN-T∂3FTt,xt,L^tA⊗In-1∂1FTt,xt,L^t+∂2FTt,xt,L^t1N⊗ft,xt,0-koL^t-FTt,xt,L^t.Clearly, solving (19) in order to obtain (17) requires invertibility of the first factor.We then defineGT:Rn×Rn×RnN→RnN,(20)GTt,x,L^=InN-T∂3FTt,x,L^A⊗In-1∂1FTt,x,L^+∂2FTt,x,L^1N⊗ft,x,0-koL^-FTt,x,L^.From (16a), (16b), (17), and (20) and denoting l^i(t) as the components of L^(t), the control law yields the following closed-loop system of ODEs:(21)x˙tL^˙t=ft,xt,Kt,xtTCL^tGTt,xt,L^tx0given,L^0=L0.The solution of (21) is an approximation of the solution of the PDE given in (7a) and (7b) and therefore only an approximate stabilization of the initial orbit is expected or rather the stabilization to an orbit close to the initial one.The ODE (21) has two types of state components, corresponding to the controlled system dynamics and to the dynamical state controller. Once L^(t) stands for a set of unmeasured state variable components, (17) can be interpreted as a state observer. Notice that this estimator introduces dynamical feedback whose state has a dimension equal to the number of points of the R-K method adopted multiplied by the dimension of the initial system to be controlled (nN).


## 3.1. Approximation of the Prediction Term: 1st Step

The first step consists in approximating the solution of (7b) by an implicit R-K ODE integration method [32], in order to estimate the prediction term, that is, the terminal value,(9)Xt,T=x+T∫01ft+Ts,Xt,Ts,0ds,xgiven.In order to estimateX(t,T) given by (9), the state transition map of the free system φ(t2,t1,x,0) is first approximated by the operator z defined by(10a)zt+T,t,x=x+T∑i=1Ncilit(10b)lit=ft+Tsi,x+T∑j=1Naijljt,0,where lj:R→Rn and aij,ci∈R are weights chosen according to the implicit method used [32]. The approximation z(t+T,t,x(t)) of X(t,Ts),  s∈[0,1], is calculated at the discretization points s=si,  i=1,…,N. N shall be chosen accordingly such that the sums in (10a) and (10b) lead to an adequate approximation of the integral in (9) for the given f(t,x,u).For simplicity, (10b) is written in the vector form (11)Lt=FTt,x,Lt,where (12)Lt=l1t⋮lNt∈RnN,∀t≥0, and FT:R×Rn×RnN→RnN is defined by (13)FTt,x,L=ft+Ts1,x+∑j=1Na1jlj,0⋮ft+TsN,x+∑j=1NaNjlj,0,∀x∈Rn,∀L∈RnN.In order to computez(t+T,t,x) through (10a), it is necessary to solve the algebraic system of (11) with unknown L(t)∈RN.Writing(14)C=c1⋯cN and closing (5) by (15)ut,xt=Kt,xtzt+T,t,xt-xt=TKt,xtCLt yield the differential algebraic equation (DAE)(16a)x˙t=ft,xt,TKt,xtCLt,x0given(16b)Lt=FTt,xt,Lt.The real-time solution of the DAE (16a) and (16b) requires the computation of L(t), its algebraic term, at each time t. We therefore introduce an observer equation in the sequel, to transform the controlled system into a system of ODEs.


## 3.2. Approximation of the Prediction Term: 2nd Step

We now approximate (16b) by solving the nN-dimensional ODE (17) whose solution L^(t) is an estimation of L(t). The initial value L^(0) is intended to be (precisely) computed off-line to provide good tracking quality for L(t).(17)ddtL^t-FTt,xt,L^t+koL^t-FTt,xt,L^t=0,L^0=L0.The scalar gainko is chosen positive in order that L^(t) of (17) tends asymptotically towards the solution L(t) of (16b) when t→+∞ and typically in such a way that the estimator dynamics are faster than the controlled system dynamics. If the evolution of L^(t) may indeed be chosen in order to fulfill (17), convergence does occur.From (17), we obtain(18)L^˙t-∂1FTt,xt,L^t+∂2FTt,xt,L^t1N⊗ft,xt,0+∂3FTt,xt,L^tTA⊗InL^˙t+koL^t-FTt,xt,L^t=0,where A=(aij) and ∂i is the partial derivative with respect to the ith variable. In is the n×n identity matrix and 1N is the column vector of dimension N with all elements equal to 1. ⊗ is the Kronecker product.By determiningL^˙(t) from (18) we get (19)L^˙t=InN-T∂3FTt,xt,L^tA⊗In-1∂1FTt,xt,L^t+∂2FTt,xt,L^t1N⊗ft,xt,0-koL^t-FTt,xt,L^t.Clearly, solving (19) in order to obtain (17) requires invertibility of the first factor.We then defineGT:Rn×Rn×RnN→RnN,(20)GTt,x,L^=InN-T∂3FTt,x,L^A⊗In-1∂1FTt,x,L^+∂2FTt,x,L^1N⊗ft,x,0-koL^-FTt,x,L^.From (16a), (16b), (17), and (20) and denoting l^i(t) as the components of L^(t), the control law yields the following closed-loop system of ODEs:(21)x˙tL^˙t=ft,xt,Kt,xtTCL^tGTt,xt,L^tx0given,L^0=L0.The solution of (21) is an approximation of the solution of the PDE given in (7a) and (7b) and therefore only an approximate stabilization of the initial orbit is expected or rather the stabilization to an orbit close to the initial one.The ODE (21) has two types of state components, corresponding to the controlled system dynamics and to the dynamical state controller. Once L^(t) stands for a set of unmeasured state variable components, (17) can be interpreted as a state observer. Notice that this estimator introduces dynamical feedback whose state has a dimension equal to the number of points of the R-K method adopted multiplied by the dimension of the initial system to be controlled (nN).


## 4. aPBC Application Issues

Two aspects of the aPBC shall be considered for application. One is the choice of an implicit R-K method to solve (10a) and (10b) and the other is the control gain tuning for stabilization. Here, we use the orthogonal collocation method with Lagrange polynomials and a constant control gain K computed through the closed-loop monodromy matrix described in the sequel.


### 4.1. Orthogonal Collocation Method with Lagrange Polynomials

The implicit R-K method given in (10a) and (10b) is a general formulation used for the integration of differential equations whose application depends on the choice of a specific implementation. Herein, using the results in [28], we apply the orthogonal collocation [30–32] as implicit R-K method.Collocation methods amount to approximate the prediction term byz(t+T,t,x(t)), where z(t+Ts,t,x(t)) is defined on the whole interval s∈[0,1] by(22)zt+Ts,t,xt=∑j=1Nwjsmjt,m1t=xt,s∈0,1,t≥0,where mj:R+→Rn,  j=1,…,N, is a column vector. The functions wj(s) are the Lagrange polynomials(23)wjs=∏i=1,i≠jNs-sisj-si,j=1,…,N,attached to the choice of the points 0=s1<s2<⋯<sN-1<sN=1.The link with the implicit R-K method is as follows.Theorem 1 (adapted from Theorem7.7 in [32]).
The collocation method (22) with Lagrange polynomials wj(s) is equivalent to the N-stage implicit Runge-Kutta method (10a) and (10b) with coefficients(24)aij=∫0siwjsds,cj=∫01wjsds,i,j=1,…,N.It is possible to choose the collocation points in order to fulfill the orthogonality relations:(25)∫011-sswiswjsds=0,i,j=2,…,N-1,i≠j.Note that, for eachj=1,…,N,  φ(t+Tsj,t,x(t),0)≈z(t+Tsj,t,x(t))=mj(t) as wj(si)=δji (with δji, the Kronecker delta symbol). A characteristic of the orthogonal collocation method is that each mj(t) is an approximation of the state that differs from lj(t) in (10a) and (10b), which are the derivatives.Equation (21) is obtained applying the substitutions provided by (22) and Theorem 1 with interpolating times sj obtained by solving (25).


### 4.2. Closed-Loop Monodromy Matrix

The aPBC was developed for the general control gainK(t,x(t)) of system (21). However, methods to choose the time and state dependent control gain have not yet been studied. Moreover, the control gain is tuned for the PBC method and applied to the aPBC, justified by the precision of the estimation of future states [28]. For simplicity, a constant control gain K is applied that depends upon the ability of computing the closed-loop monodromy matrix of x∗(t). The computation of this matrix requires the integration of the closed-loop system and its variational equation (26) along a trajectory in the vicinity of x∗(t) (linearised dynamics around the periodic orbit) [33, Appendix B]. To integrate this trajectory, the initial condition is chosen close to x∗(0). Integrating (26) over a period yields the corresponding closed-loop monodromy matrix Ψ(t)=Φ(t+T,t).(26)dΦt,0dt=∇xft,x,Kφt+T,t,x,0-xx=x∗tΦt,0,t∈0,T,Φ0,0=In,where(27)∇xft,x,Kφt+T,t,x,0-x=∇xft,x,u+K∇uft,x,u∇xφt+T,t,x,0-In.∇xφ(t+T,t,x,0) is the free system monodromy matrix.Using (26) and (27), we compute the closed-loop monodromy matrix of x∗(t) given gain K. The Floquet multipliers are computed to measure the local stability of the controlled orbit for the chosen K.In practice, we fixK, compute the monodromy matrix by integrating (26) with an explicit R-K method, and φ(t+T,t,x∗(t),0) is computed by integrating the free system over a period T at each step of the integration of (26). After that, we obtain the corresponding Floquet multipliers of the closed-loop system.


## 4.1. Orthogonal Collocation Method with Lagrange Polynomials

The implicit R-K method given in (10a) and (10b) is a general formulation used for the integration of differential equations whose application depends on the choice of a specific implementation. Herein, using the results in [28], we apply the orthogonal collocation [30–32] as implicit R-K method.Collocation methods amount to approximate the prediction term byz(t+T,t,x(t)), where z(t+Ts,t,x(t)) is defined on the whole interval s∈[0,1] by(22)zt+Ts,t,xt=∑j=1Nwjsmjt,m1t=xt,s∈0,1,t≥0,where mj:R+→Rn,  j=1,…,N, is a column vector. The functions wj(s) are the Lagrange polynomials(23)wjs=∏i=1,i≠jNs-sisj-si,j=1,…,N,attached to the choice of the points 0=s1<s2<⋯<sN-1<sN=1.The link with the implicit R-K method is as follows.Theorem 1 (adapted from Theorem7.7 in [32]).
The collocation method (22) with Lagrange polynomials wj(s) is equivalent to the N-stage implicit Runge-Kutta method (10a) and (10b) with coefficients(24)aij=∫0siwjsds,cj=∫01wjsds,i,j=1,…,N.It is possible to choose the collocation points in order to fulfill the orthogonality relations:(25)∫011-sswiswjsds=0,i,j=2,…,N-1,i≠j.Note that, for eachj=1,…,N,  φ(t+Tsj,t,x(t),0)≈z(t+Tsj,t,x(t))=mj(t) as wj(si)=δji (with δji, the Kronecker delta symbol). A characteristic of the orthogonal collocation method is that each mj(t) is an approximation of the state that differs from lj(t) in (10a) and (10b), which are the derivatives.Equation (21) is obtained applying the substitutions provided by (22) and Theorem 1 with interpolating times sj obtained by solving (25).


## 4.2. Closed-Loop Monodromy Matrix

The aPBC was developed for the general control gainK(t,x(t)) of system (21). However, methods to choose the time and state dependent control gain have not yet been studied. Moreover, the control gain is tuned for the PBC method and applied to the aPBC, justified by the precision of the estimation of future states [28]. For simplicity, a constant control gain K is applied that depends upon the ability of computing the closed-loop monodromy matrix of x∗(t). The computation of this matrix requires the integration of the closed-loop system and its variational equation (26) along a trajectory in the vicinity of x∗(t) (linearised dynamics around the periodic orbit) [33, Appendix B]. To integrate this trajectory, the initial condition is chosen close to x∗(0). Integrating (26) over a period yields the corresponding closed-loop monodromy matrix Ψ(t)=Φ(t+T,t).(26)dΦt,0dt=∇xft,x,Kφt+T,t,x,0-xx=x∗tΦt,0,t∈0,T,Φ0,0=In,where(27)∇xft,x,Kφt+T,t,x,0-x=∇xft,x,u+K∇uft,x,u∇xφt+T,t,x,0-In.∇xφ(t+T,t,x,0) is the free system monodromy matrix.Using (26) and (27), we compute the closed-loop monodromy matrix of x∗(t) given gain K. The Floquet multipliers are computed to measure the local stability of the controlled orbit for the chosen K.In practice, we fixK, compute the monodromy matrix by integrating (26) with an explicit R-K method, and φ(t+T,t,x∗(t),0) is computed by integrating the free system over a period T at each step of the integration of (26). After that, we obtain the corresponding Floquet multipliers of the closed-loop system.


## 5. Case Studies

Two continuous-time systems are used as case study for the aPBC—the forced van der Pol oscillator and Rössler system. In this section we provide a brief analysis about them evidencing their chaotic behavior through bifurcation diagrams and chaotic attractors.


### 5.1. Forced van der Pol Oscillator

The forced van der Pol (vdP) oscillator is described by the nonlinear state space model (28):(28)x˙t=x2tasin⁡ωt-ηx12t-1x2t-x1t,where x:R+→R2 and f:R+×R2→R2. In the present work we set ω=0.45 and η=1 and provide bifurcation diagrams by attractor Poincaré map (Figure 1(a)) and largest nonzero Lyapunov exponent (λmax⁡, Figure 1(b)) for a=0.9850.989 evidencing the route to chaos. Poincaré section is defined as usual for nonautonomous systems using time, and in this case (29)SvdP=svdP=x1,x2∈R2:tmod⁡T=0,T=2πω, where (Amod⁡B) is the modulo operation or the remainder of the integer division of A by B.Figure 1
Bifurcation diagram of forced van der Pol oscillator forω=0.45 and η=1 using (a) Poincaré map and (b) largest nonzero Lyapunov exponent.
(a)
(b)Form Figure1 it is possible to note that increasing a we have a series of period-doubling bifurcations leading to two chaotic attractors that are merged by a merging crises bifurcation. Lyapunov exponents are computed for each a evidencing the bifurcations and chaotic attractors. This work follows [29] and adopts a=0.988 with chaotic attractor in continuous-time and Poincaré map shown in Figure 2.Figure 2
Forced van der Pol chaotic attractor in continuous-time (CA) and Poincaré map (CA-PM) forω=0.45, η=1, and a=0.988.


### 5.2. Rössler System

The Rössler system is described by the nonlinear state space model(30)x˙t=-x2t-x3tx1t+β1x2tβ2+x1tx3t-β3x3t,where x:R+→R3 and f:R+×R3→R3. Here, β1=β2=0.2 and the route to chaos through period-doubling bifurcations for β3=[25] is presented by attractors Poincaré map (Figure 3(a)) and largest Lyapunov exponent (λmax⁡, Figure 3(b)) bifurcation diagrams. This work follows [34] and adopts β3=4.5 with chaotic attractor in continuous-time and Poincaré map shown in Figure 4. Poincaré section is defined for the Rössler system in (31)SRs=sRs=x1,x2,x3∈R3:x1=0,x2≥0.Figure 3
Bifurcation diagram of Rössler system forβ1=β2=0.2 using (a) Poincaré map and (b) largest Lyapunov exponent.
(a)
(b)Figure 4
Rössler system chaotic attractor in continuous-time (CA) and Poincaré map (CA-PM) forβ1=β2=0.2 and β3=4.5.


## 5.1. Forced van der Pol Oscillator

The forced van der Pol (vdP) oscillator is described by the nonlinear state space model (28):(28)x˙t=x2tasin⁡ωt-ηx12t-1x2t-x1t,where x:R+→R2 and f:R+×R2→R2. In the present work we set ω=0.45 and η=1 and provide bifurcation diagrams by attractor Poincaré map (Figure 1(a)) and largest nonzero Lyapunov exponent (λmax⁡, Figure 1(b)) for a=0.9850.989 evidencing the route to chaos. Poincaré section is defined as usual for nonautonomous systems using time, and in this case (29)SvdP=svdP=x1,x2∈R2:tmod⁡T=0,T=2πω, where (Amod⁡B) is the modulo operation or the remainder of the integer division of A by B.Figure 1
Bifurcation diagram of forced van der Pol oscillator forω=0.45 and η=1 using (a) Poincaré map and (b) largest nonzero Lyapunov exponent.
(a)
(b)Form Figure1 it is possible to note that increasing a we have a series of period-doubling bifurcations leading to two chaotic attractors that are merged by a merging crises bifurcation. Lyapunov exponents are computed for each a evidencing the bifurcations and chaotic attractors. This work follows [29] and adopts a=0.988 with chaotic attractor in continuous-time and Poincaré map shown in Figure 2.Figure 2
Forced van der Pol chaotic attractor in continuous-time (CA) and Poincaré map (CA-PM) forω=0.45, η=1, and a=0.988.


## 5.2. Rössler System

The Rössler system is described by the nonlinear state space model(30)x˙t=-x2t-x3tx1t+β1x2tβ2+x1tx3t-β3x3t,where x:R+→R3 and f:R+×R3→R3. Here, β1=β2=0.2 and the route to chaos through period-doubling bifurcations for β3=[25] is presented by attractors Poincaré map (Figure 3(a)) and largest Lyapunov exponent (λmax⁡, Figure 3(b)) bifurcation diagrams. This work follows [34] and adopts β3=4.5 with chaotic attractor in continuous-time and Poincaré map shown in Figure 4. Poincaré section is defined for the Rössler system in (31)SRs=sRs=x1,x2,x3∈R3:x1=0,x2≥0.Figure 3
Bifurcation diagram of Rössler system forβ1=β2=0.2 using (a) Poincaré map and (b) largest Lyapunov exponent.
(a)
(b)Figure 4
Rössler system chaotic attractor in continuous-time (CA) and Poincaré map (CA-PM) forβ1=β2=0.2 and β3=4.5.


## 6. Numerical Results

Numerical results are presented providing illustrative examples about the aPBC characteristics. We use the forced van der Pol oscillator as case study to show an UPO that cannot be stabilized by the DFC method due to the odd-number limitation and is stabilized using the aPBC. Using the same dynamical system, the robustness of the aPBC is evaluated for low-precision future states estimation. The aPBC is applied to the autonomous Rössler system and its robustness is also evaluated. For all the cases, the valuesN, ko, and K that characterize the aPBC are tuned prior to application.


### 6.1. Applying aPBC to the Forced van der Pol Oscillator

In this section, we provide numerical results on the approximated Prediction-Based Control (aPBC) applied to the forced van der Pol (vdP) oscillator according to (32), where f(t,x(t)) is defined in (28).(32)x˙t=ft,xt+ut,xt.The state transition mapφ(t+Tsj,t,x(t),0) in (4) is approximated by an implicit R-K method using the operator z(t+Tsj,t,x(t)) in the sj,  j=1,2,…,N, discretization points. The implicit R-K method is solved using estimator (17) resulting in the aPBC (21).Herein, we apply the aPBC using the orthogonal collocation method (see Section4.1), obtaining z(t+Tsj,t,x(t))=mj(t), and estimator (17) that leads to m^j(t)≈mj(t). The control signal of the closed-loop control method (21) applied to system (32) is(33)ut,xt=Km^Nt-xtwith the constant control gain(34)K=k000,k∈R.u:R+×R2→R2 is defined with the same dimension of the free system, but only the first component of the control signal vector is different from zero. This choice of matrix K makes only the first state variable influence on feedback and simplifies its tuning; however, better results can be obtained for a full matrix K.The control signal is designed by approximating the feedback term, leading to the closed-loop system (21). The solution of the latter depends on the choice of three parameters: the control gain k, the observer gain ko, and the number of discretization points for the orthogonal collocation method N. The values ko and N are directly related to the estimation quality and these two parameters will be tuned first.


#### 6.1.1. Tuning the aPBC Parametersko and N

In order to evaluate the future state estimation, whose characteristics are related to the parametersko and N, we first set k=0. The system is integrated for t∈[0,T] and the points m^j(t+jT/N), j=0,…,N, are collected for different values of ko and N. One then assesses the convergence by comparison with the solution obtained with N=152 and ko=50, by computing (35)εN,ko=m^j,i152,50·-m^j,iN,ko·L10,Tm^j,i152,50·L10,T.The superscript on m^j,i(N,ko)(t) indicates the value of N and ko used in the estimation and the subscript i refers to the state variable of the n-dimensional vector m^j(t).Figure5(a) shows ε(N,50) for different values of N using initial condition x(0)=-0.750.75T,  T=2π/0.45, and state variable i=2. Note that the collocation parameters si and ωj are previously computed for each N tested. The computation of ωj for large N requires intense off-line computational burden and reduces the number of values N tested. These values are computed once for all and are independent of f(t,x(t)). The values of N<72 lead to integration instability and are not computed.Figure 5
Relative error of the free system response estimated future value using different (a)N and (b) ko for the forced vdP oscillator.
(a)
(b)Repeating the same process, we evaluateε(152,ko) for different values of ko; see Figure 5(b). Increasing N and ko results in larger computation effort and trade-off between estimation quality and computation effort should be considered.Due to the error level shown in Figure5 for N=102 and ko=10, we adopt these values in the sequel. Figure 6 shows x2(t) and m^j,2(t-T) (beware of the time shift here). In this case, it is expected that x(t)=φ(t,t-T,x(t-T),0)≈z(t,t-T,x(t-T))≈m^N(102,10)(t-T) and the proximity of both time-series in Figure 6 allows the use of mj(t) instead of m^j(102,10)(t). Once the computational cost of the aPBC comes from processing the set of nN ODEs used for prediction with reasonable error tolerance, m^j(102,10)(t) is also a reference for a good trade-off between future state estimation and computation cost.Figure 6
Current time and shifted future time trajectories of the forced vdP oscillator fork=0, ko=10, and N=102.


#### 6.1.2. Target Unstable Periodic Orbits

Two UPOs are chosen as target UPO (TUPO) to be stabilized by the aPBC. TUPO1 was used in [29] for the application of an optimal DFC and in [28] for comparison between aPBC and optimal DFC. Herein, it is used to exemplify the effect of low-precision future state estimation through a small number of collocation points N. TUPO2 is used to exemplify aPBC stabilizing an orbit subjected to the odd-number limitation when applying DFC.Figure7 shows the state space of the target periodic orbits, in continuous-time (line) and Poincaré map (∗), embedded in the chaotic attractor. One point of the orbits with adequate precision for open-loop simulations and their Floquet multipliers are presented in Table 1. Period-1 (where period-1 means one cycle of the forced term or a one-point Poincaré map for the nonautonomous system) orbits of the forced van der Pol oscillator with ω=0.45 leads to T=2π/0.45.Table 1
Floquet multipliers and one Poincaré map point fora=0.988 of the two target periodic orbits of the vdP oscillator.
TUPO1
TUPO2
x
∗
(
t
0
)
[
0.1588445
-
0.1106056
]
T
[
0.0771541
1.6893651
]
T
μ
[−1.872 −0.048]
[
4.627
0.147
]Figure 7
Target UPOs (TUPO1 and TUPO2) embedded in the chaotic attractor (CA) of the forced vdP oscillator. The points for the target UPOs and the chaotic attractor in the Poincaré map (PM) are represented by asterisk and gray bullet, respectively.
(a)
(b)


#### 6.1.3. Stabilizing TUPO1 through aPBC for Low-Precision Future State Estimation

In Section6.1.1, the values of N=102 and ko=10 that result in a good estimation of future states according to value ε(N,ko) have been obtained. Results for these N and ko and a comparison with DFC are shown in [28]. In the present section, the aPBC is applied to stabilize TUPO1 using N<102 exemplifying the robustness of the method to the use of low-precision future state estimates, aiming at possible experimental applications.Before applying aPBC to TUPO1, a control gaink that leads to closed-loop stabilization shall be found.Tuning ofk.We use (26), with x∗(0)=0.1588445-0.1106056T, to obtain a value k that stabilizes the orbit of (7a) and (7b). Note that, once the estimation with N=102 and ko=10 is precise, a stabilizing k for (7a) and (7b) yields stability for a slightly perturbed cycle of (32) with the control signal (33) [28].Figure8 shows the magnitude of the two Floquet multipliers |μ| for different values of k. Stability is achieved when |μi|≤1,i=1,2. We choose k=0.25, which results in |μ|max⁡≈0.01213. All the Floquet multipliers in this case are real numbers.Figure 8
Magnitude of the two Floquet multipliers for different values ofk for the TUPO1 controlled by the aPBC.Applying the aPBC for a LowN.Figure 5(a) shows an increase of ε(N,50) when decreasing N for k=0. This value refers to the free system and now we are then interested in characteristics of the closed-loop system, particularly the value of steady state error and control effort as a function of N.Steady state error indexd is the distance (Euclidean norm) between predicted state m^N(t) and current state x(t) for an entire cycle in the close vicinity of target UPO. Steady state control effort index υ is the measure of total external signal (norm-1) necessary to keep the trajectory in close vicinity around one cycle of the target periodic orbit.(36)d=limt→+∞⁡m^Nτ-xτL2t,t+T.υ=limt→+∞⁡uτL1t,t+T.Figure9(a) shows d when stabilizing the TUPO1 by applying the aPBC for different values of N with k=0.25 and ko=10. Figure 9(b) shows the same for υ.Figure 9
Indexes (a)d and (b) υ when applying aPBC to TUPO1 for different N with k=0.25 and ko=10. Bottom dashed line are the values of (a) d and (b) υ computed with N=152 and ko=50 used as reference.
(a)
(b)It can be observed in Figure9 that even for a reduction in the number of collocation points, stability of (an orbit close to) TUPO1 is achieved for N≥67. TUPO1 stabilized with the aPBC for N=82 and N=67 and x(0)=-1.51.0T are shown in Figures 10 and 11, respectively. Figures 10(e) and 11(e) show the initial value for the predicted states mj(0), j=1,2,…,N, and φ(T,0,x(0),0). It is observed that mj(0) is relatively close to φ(T,0,x(0),0).Figure 10
aPBC applied withN=82, k=0.25, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c) and (d)) time-series of the actual and predicted state variables; (e) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)Figure 11
aPBC applied withN=67, k=0.25, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c) and (d)) time-series of the actual and predicted state variables; (e) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)The high values observed in Figure9 for N=3 are explained in Figure 12. It exemplifies a situation wherein a stable orbit different from TUPO1 is achieved due to the reduced N. High amplitude periodic control signal in Figure 12(b) indicates that this stabilized orbit is not an orbit of the original free system. Figure 12(e) shows that N=3 leads to mj(0) almost unrelated to φ(T,0,x(0),0) and consequently u2(t)=k(m^N,2(t)-x2(t)) is a high amplitude control signal. According to Figure 12(a), a trajectory with initial condition x(0) on TUPO1 (Poincaré point) diverges from it showing that it is not stable in this case.Figure 12
aPBC applied withN=3, k=0.25, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c) and (d)) time-series of the actual and predicted state variables; (e) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)It is noticed from the time-series of trajectories in Figures10(c), 10(d), 11(c), and 11(d) that in certain torsions of the trajectory of the system controlled by the aPBC there is an oscillation on m^N(t). This oscillation is generated because small N leads to instability on the predicted state and it reflects in time-series of control signals in Figures 10(b) and 11(b) and on the values of υ and d in Figure 9.Another consequence on reducingN is the loss of stability of m^N(t), exemplified in Figure 13 and occurring even after changing the error tolerances of the explicit R-K integrator (the simulations present in this work were performed using the MATLAB® integration routine ode45 and, in this specific case, the routines ode113 and ode15s were also tried, all of them resulting in interruption of the numerical integration). One of the reasons for this loss of stability is that N is the parameter that rules the implicit R-K integration error and Figure 13(d) shows that mN(0) slightly diverges from φ(T,0,x(0),0) at time T and, contrasting to example of Figure 12, the trajectories of the new system do not converge to a stable periodic orbit.Figure 13
Time-series of a trajectory when applying the aPBC withN=62 and (a) k=0 and ko=10; (b) k=0.25 and ko=10; (c) k=0.25 and ko=1000; (d) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)


#### 6.1.4. Stabilizing TUPO2 through aPBC: An Odd-Number Limitation Example

One of the main purposes when developing PBC for discrete-time systems is the odd-number limitation subjected by the DFC [14, 21]. Herein, we show a case where the continuous-time DFC with a constant control gain is subjected to the odd-number limitation and aPBC is not. TUPO2 has an odd number of real Floquet multipliers larger than +1 (see Table 1) and is used for a numerical example.Applying the DFC.The DFC with control signal (37)ut,xt=Kxt-T-xt,K∈R2×2, is applied to system (32) with the goal of stabilizing TUPO2. The procedure of finding K for DFC follows [29].A scalar control gaink is used with K defined in (34) resulting in the root locus chart of Figure 14. This chart contains the free system orbit Floquet multipliers μ marked with ×: in black solid line, the fourth largest magnitude (the closed-loop system is infinite-dimensional when applying DFC; see [29] for details) closed-loop Floquet multipliers of the orbit for positive k and, in dark gray, the largest closed-loop Floquet Multiplier of the orbit for negative k. Arrows indicate increasing |k|. It is possible to conclude that for any k there is a Floquet multiplier with magnitude larger than 1 and the orbit is not stabilizable. In addition, all the Floquet multipliers of the orbit tend to +1 when increasing DFC positive scalar gain k. The same chart was done for the other three possible scalars k with equivalent (omitted) results.Figure 14
Root locus chart for the Floquet multipliers of the TUPO2 using the DFC. In black, the four largest magnitude Floquet multipliers for positive control gaink. In dark gray, the largest magnitude Floquet multiplier for negative control gain k. The Floquet multipliers of the free system orbit are represented by (×) and the unit circle by dashed line in black.Following [29], a matrix gain was tried with values obtained by an optimization process using the MATLAB routinesfminsearch andfmincon [35–38]. This optimization aims at finding a minimum for the largest magnitude Floquet multiplier through changing matrix K. A scan process varying all the values of the matrix gain was also tried. The final trial used a matrix gain composed of a rotational matrix (formed by sines and cosines) multiplied by a real scalar gain equivalent to the gain used in [18–20], varying both the rotation angle and the gain. After performing the three methodologies, a matrix control gain that stabilizes TUPO2 was not found.Applying aPBC and Stabilizing TUPO2.Figure 15 shows the magnitude of the Floquet multipliers as a function of the control gain k for the TUPO2 by applying the aPBC. It was noticed that the aPBC stabilizes orbits with an odd number of real Floquet multipliers larger than +1. One detail is relevant here; the stabilization is achieved only for negative values of the gain k and positive values of k make the largest magnitude Floquet multiplier increase.Figure 15
Magnitude of the Floquet multipliers for different values ofk for the TUPO2 controlled by the aPBC.TUPO2 stabilized by the aPBC withk=-0.125 and initial condition x(0)=0.11.8T is shown in state space in Figure 16(a). Time-series of the control signal is presented in Figure 16(b) with u2(t)→0 as t→+∞ indicating stabilization.Figure 16
TUPO2 of forced vdP oscillator stabilized using aPBC withk=-0.125 and x0=[0.11.8]T. (a) State space trajectory and (b) time-series of the control effort.
(a)
(b)


### 6.2. Applying aPBC to the Rössler System

In this section, we provide numerical results on approximated Prediction-Based Control (aPBC) for the Rössler system using control scheme (21) applied to system (30) according to (32). The goals are (i) to show that the method is applicable to autonomous systems without any modification and (ii) to compare the effect of different numbers of collocation points N on stabilization.aPBC is applied following the same procedure shown in Section6.1 with control signal u(t,x(t)) defined in (33). Herein, we use the constant control gain(38)K=00000000k,k∈R.This choice of matrix K makes only the third state variable influence on feedback and simplifies its tuning.


#### 6.2.1. Tuning the aPBC Parametersko and N

The procedure shown in Section6.1.1 is repeated here to characterize the performance of the estimator for different values of the parameters ko and N for k=0. We used initial condition x(0)=0-60.0375T to compute ε(N,ko) in (35). The results are shown in Figure 17 for different values of (a) N with ko=50 (ε(N,50)) and (b) ko with N=152 (ε(152,ko)).Figure 17
Relative error of the free system response estimated future value using different (a)N and (b) ko for the Rössler system.
(a)
(b)In this caseN<32 values lead to numerical integration instability, except when N=3, and are not computed. Herein we notice that less collocation points are necessary to obtain estimation quality equivalent to the one obtained for the vdP oscillator (see Figure 5). It can be seen that for the vdP there is ε(102,50)≈4.4×10-4, while for the Rössler system there is ε(52,50)≈4.4×10-4. Values ε(152,ko) for large ko of the chosen trajectory of the Rössler system also stabilize to a value much smaller than for the vdP oscillator.Figure18 shows x2(t) and m^N,2(t-T) (beware of the time shift here) for the Rössler system, which can be compared with Figure 6. Notice that T=5.8439698764 (the period of the TUPO4) for the Rössler system while T=2π/0.45 for vdP oscillator, but the reduction on ε(152,ko) is much more related to the behavior of the trajectories than the period itself, because collocation points si are computed for a normalized time between 0 and 1.Figure 18
Current time and shifted future time trajectories of the Rössler system fork=0, ko=10, and N=102.


#### 6.2.2. Target Unstable Periodic Orbit

Asenjo et al. [34] present one UPO of the chaotic attractor for β1=β2=0.2 and β3=4.5. We apply the Newton-Raphson method with Poincaré section (31) to refine the position and period of the orbit and call it TUPO4 (Figure 19). One point of TUPO4 and its period T with adequate precision for open-loop simulations and Floquet multipliers μ are presented in Table 2.Table 2
Poincaré map pointx(t0), period T, and Floquet multipliers μ of a periodic orbit (TUPO4) of the Rössler system.
TUPO4
x
(
t
0
)
0
5.0206004746
1.6806432235
T
T
5.8439698764
μ
-
1.918
1
-
1.866
×
1
0
-
11Figure 19
TUPO4 of Rössler system in continuous-time (line) and Poincaré map point (∗) for β1=β2=0.2 and β3=4.5.


#### 6.2.3. Stabilizing TUPO4 through aPBC

This section is dedicated to the stabilization of TUPO4 using aPBC for different valuesN. From Figure 17, values N=102 and ko=10 are indicated for an accurate estimation of future states and aPBC will be initially applied for these values. Results for low N are also presented, evidencing the robustness of the controller for low-precision future states estimation in autonomous systems.The control gaink of (38) shall be tuned for TUPO4 prior to the application.Tuning ofk.Figure 20 shows the magnitude of the three Floquet multipliers |μ| for different values of k according to the procedure presented in Section 6.1.3. The value μ1=1 is characteristic of autonomous system trajectories and stability is achieved when, for the other Floquet multipliers, |μi|≤1,i=2,3. Choosing k=1.85 results in |μ2|≈2.994×10-3 (and |μ3|≈2.167×10-14). All Floquet multipliers in this case are real numbers.Figure 20
Magnitude of the Floquet multipliers for different values ofk for the TUPO4 controlled by the aPBC.Applying the aPBC for Different Values ofN.In this section, we evaluate the characteristics of aPBC for different values N applied to the stabilization of TUPO4. Herein, we use N≤102,  ko=10, and k=1.85.Figure17(a) shows an increase of ε(N,50) observed when decreasing N for k=0. Figure 21 shows d and υ for different values of N with k=1.85 and ko=10. Values of 3<N<32 were not computed because of instability, the same as observed in Section 6.1.3.Figure 21
Indexes (a)d and (b) υ when applying aPBC to TUPO4 for different N with k=1.85 and ko=10. The bottom dashed lines are the values of (a) d and (b) υ computed with N=152 and ko=50 used as reference.
(a)
(b)Herein, we observe thatN≥102 leads to a precise estimation once there is a convergence of υ to the value computed using N=152 and ko=50 (dashed line). Comparing Figure 21 with Figure 9, a great reduction on υ and d is observed for the control method applied to the Rössler system, which is in accordance with what was observed for ε(N,50). In fact, N between 47 and 52 leads to a value of d equivalent to the value obtained for N=102 for the vdP oscillator. Considering υ, this reduction is more significant, with N=37 leading to a value of υ close to the value obtained for N=102 for the vdP. This result is important for practical aspects because it shows that, for certain systems, the augmented state of the closed-loop system controlled with the aPBC can be reduced significantly (in comparison with the initial results shown for the vdP oscillator).Figure21 shows that stability of TUPO4 is achieved for N≥32. The TUPO4 stabilized by the aPBC for N=102, N=52, and N=32 is shown in Figures 22, 23, and 24, respectively, using x(0)=0-60.0375T. Deviation of m^N,2(t) from x2(t) is observed in Figure 24(d), which reflects the numerical integration instabilities for N<32. We verify that stabilization is not achieved for N=3, which leads to the high values of υ and d observed in Figure 21.Figure 22
aPBC applied to the Rössler system withN=102, k=1.85, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c), (d), and (e)) time-series of the actual and predicted state variables; (f) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)
(f)Figure 23
aPBC applied to the Rössler system withN=52, k=1.85, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c), (d), and (e)) time-series of the actual and predicted state variables; (f) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)
(f)Figure 24
aPBC applied to the Rössler system withN=32, k=1.85, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c), (d), and (e)) time-series of the actual and predicted state variables; (f) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)
(f)


## 6.1. Applying aPBC to the Forced van der Pol Oscillator

In this section, we provide numerical results on the approximated Prediction-Based Control (aPBC) applied to the forced van der Pol (vdP) oscillator according to (32), where f(t,x(t)) is defined in (28).(32)x˙t=ft,xt+ut,xt.The state transition mapφ(t+Tsj,t,x(t),0) in (4) is approximated by an implicit R-K method using the operator z(t+Tsj,t,x(t)) in the sj,  j=1,2,…,N, discretization points. The implicit R-K method is solved using estimator (17) resulting in the aPBC (21).Herein, we apply the aPBC using the orthogonal collocation method (see Section4.1), obtaining z(t+Tsj,t,x(t))=mj(t), and estimator (17) that leads to m^j(t)≈mj(t). The control signal of the closed-loop control method (21) applied to system (32) is(33)ut,xt=Km^Nt-xtwith the constant control gain(34)K=k000,k∈R.u:R+×R2→R2 is defined with the same dimension of the free system, but only the first component of the control signal vector is different from zero. This choice of matrix K makes only the first state variable influence on feedback and simplifies its tuning; however, better results can be obtained for a full matrix K.The control signal is designed by approximating the feedback term, leading to the closed-loop system (21). The solution of the latter depends on the choice of three parameters: the control gain k, the observer gain ko, and the number of discretization points for the orthogonal collocation method N. The values ko and N are directly related to the estimation quality and these two parameters will be tuned first.


### 6.1.1. Tuning the aPBC Parametersko and N

In order to evaluate the future state estimation, whose characteristics are related to the parametersko and N, we first set k=0. The system is integrated for t∈[0,T] and the points m^j(t+jT/N), j=0,…,N, are collected for different values of ko and N. One then assesses the convergence by comparison with the solution obtained with N=152 and ko=50, by computing (35)εN,ko=m^j,i152,50·-m^j,iN,ko·L10,Tm^j,i152,50·L10,T.The superscript on m^j,i(N,ko)(t) indicates the value of N and ko used in the estimation and the subscript i refers to the state variable of the n-dimensional vector m^j(t).Figure5(a) shows ε(N,50) for different values of N using initial condition x(0)=-0.750.75T,  T=2π/0.45, and state variable i=2. Note that the collocation parameters si and ωj are previously computed for each N tested. The computation of ωj for large N requires intense off-line computational burden and reduces the number of values N tested. These values are computed once for all and are independent of f(t,x(t)). The values of N<72 lead to integration instability and are not computed.Figure 5
Relative error of the free system response estimated future value using different (a)N and (b) ko for the forced vdP oscillator.
(a)
(b)Repeating the same process, we evaluateε(152,ko) for different values of ko; see Figure 5(b). Increasing N and ko results in larger computation effort and trade-off between estimation quality and computation effort should be considered.Due to the error level shown in Figure5 for N=102 and ko=10, we adopt these values in the sequel. Figure 6 shows x2(t) and m^j,2(t-T) (beware of the time shift here). In this case, it is expected that x(t)=φ(t,t-T,x(t-T),0)≈z(t,t-T,x(t-T))≈m^N(102,10)(t-T) and the proximity of both time-series in Figure 6 allows the use of mj(t) instead of m^j(102,10)(t). Once the computational cost of the aPBC comes from processing the set of nN ODEs used for prediction with reasonable error tolerance, m^j(102,10)(t) is also a reference for a good trade-off between future state estimation and computation cost.Figure 6
Current time and shifted future time trajectories of the forced vdP oscillator fork=0, ko=10, and N=102.


### 6.1.2. Target Unstable Periodic Orbits

Two UPOs are chosen as target UPO (TUPO) to be stabilized by the aPBC. TUPO1 was used in [29] for the application of an optimal DFC and in [28] for comparison between aPBC and optimal DFC. Herein, it is used to exemplify the effect of low-precision future state estimation through a small number of collocation points N. TUPO2 is used to exemplify aPBC stabilizing an orbit subjected to the odd-number limitation when applying DFC.Figure7 shows the state space of the target periodic orbits, in continuous-time (line) and Poincaré map (∗), embedded in the chaotic attractor. One point of the orbits with adequate precision for open-loop simulations and their Floquet multipliers are presented in Table 1. Period-1 (where period-1 means one cycle of the forced term or a one-point Poincaré map for the nonautonomous system) orbits of the forced van der Pol oscillator with ω=0.45 leads to T=2π/0.45.Table 1
Floquet multipliers and one Poincaré map point fora=0.988 of the two target periodic orbits of the vdP oscillator.
TUPO1
TUPO2
x
∗
(
t
0
)
[
0.1588445
-
0.1106056
]
T
[
0.0771541
1.6893651
]
T
μ
[−1.872 −0.048]
[
4.627
0.147
]Figure 7
Target UPOs (TUPO1 and TUPO2) embedded in the chaotic attractor (CA) of the forced vdP oscillator. The points for the target UPOs and the chaotic attractor in the Poincaré map (PM) are represented by asterisk and gray bullet, respectively.
(a)
(b)


### 6.1.3. Stabilizing TUPO1 through aPBC for Low-Precision Future State Estimation

In Section6.1.1, the values of N=102 and ko=10 that result in a good estimation of future states according to value ε(N,ko) have been obtained. Results for these N and ko and a comparison with DFC are shown in [28]. In the present section, the aPBC is applied to stabilize TUPO1 using N<102 exemplifying the robustness of the method to the use of low-precision future state estimates, aiming at possible experimental applications.Before applying aPBC to TUPO1, a control gaink that leads to closed-loop stabilization shall be found.Tuning ofk.We use (26), with x∗(0)=0.1588445-0.1106056T, to obtain a value k that stabilizes the orbit of (7a) and (7b). Note that, once the estimation with N=102 and ko=10 is precise, a stabilizing k for (7a) and (7b) yields stability for a slightly perturbed cycle of (32) with the control signal (33) [28].Figure8 shows the magnitude of the two Floquet multipliers |μ| for different values of k. Stability is achieved when |μi|≤1,i=1,2. We choose k=0.25, which results in |μ|max⁡≈0.01213. All the Floquet multipliers in this case are real numbers.Figure 8
Magnitude of the two Floquet multipliers for different values ofk for the TUPO1 controlled by the aPBC.Applying the aPBC for a LowN.Figure 5(a) shows an increase of ε(N,50) when decreasing N for k=0. This value refers to the free system and now we are then interested in characteristics of the closed-loop system, particularly the value of steady state error and control effort as a function of N.Steady state error indexd is the distance (Euclidean norm) between predicted state m^N(t) and current state x(t) for an entire cycle in the close vicinity of target UPO. Steady state control effort index υ is the measure of total external signal (norm-1) necessary to keep the trajectory in close vicinity around one cycle of the target periodic orbit.(36)d=limt→+∞⁡m^Nτ-xτL2t,t+T.υ=limt→+∞⁡uτL1t,t+T.Figure9(a) shows d when stabilizing the TUPO1 by applying the aPBC for different values of N with k=0.25 and ko=10. Figure 9(b) shows the same for υ.Figure 9
Indexes (a)d and (b) υ when applying aPBC to TUPO1 for different N with k=0.25 and ko=10. Bottom dashed line are the values of (a) d and (b) υ computed with N=152 and ko=50 used as reference.
(a)
(b)It can be observed in Figure9 that even for a reduction in the number of collocation points, stability of (an orbit close to) TUPO1 is achieved for N≥67. TUPO1 stabilized with the aPBC for N=82 and N=67 and x(0)=-1.51.0T are shown in Figures 10 and 11, respectively. Figures 10(e) and 11(e) show the initial value for the predicted states mj(0), j=1,2,…,N, and φ(T,0,x(0),0). It is observed that mj(0) is relatively close to φ(T,0,x(0),0).Figure 10
aPBC applied withN=82, k=0.25, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c) and (d)) time-series of the actual and predicted state variables; (e) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)Figure 11
aPBC applied withN=67, k=0.25, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c) and (d)) time-series of the actual and predicted state variables; (e) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)The high values observed in Figure9 for N=3 are explained in Figure 12. It exemplifies a situation wherein a stable orbit different from TUPO1 is achieved due to the reduced N. High amplitude periodic control signal in Figure 12(b) indicates that this stabilized orbit is not an orbit of the original free system. Figure 12(e) shows that N=3 leads to mj(0) almost unrelated to φ(T,0,x(0),0) and consequently u2(t)=k(m^N,2(t)-x2(t)) is a high amplitude control signal. According to Figure 12(a), a trajectory with initial condition x(0) on TUPO1 (Poincaré point) diverges from it showing that it is not stable in this case.Figure 12
aPBC applied withN=3, k=0.25, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c) and (d)) time-series of the actual and predicted state variables; (e) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)It is noticed from the time-series of trajectories in Figures10(c), 10(d), 11(c), and 11(d) that in certain torsions of the trajectory of the system controlled by the aPBC there is an oscillation on m^N(t). This oscillation is generated because small N leads to instability on the predicted state and it reflects in time-series of control signals in Figures 10(b) and 11(b) and on the values of υ and d in Figure 9.Another consequence on reducingN is the loss of stability of m^N(t), exemplified in Figure 13 and occurring even after changing the error tolerances of the explicit R-K integrator (the simulations present in this work were performed using the MATLAB® integration routine ode45 and, in this specific case, the routines ode113 and ode15s were also tried, all of them resulting in interruption of the numerical integration). One of the reasons for this loss of stability is that N is the parameter that rules the implicit R-K integration error and Figure 13(d) shows that mN(0) slightly diverges from φ(T,0,x(0),0) at time T and, contrasting to example of Figure 12, the trajectories of the new system do not converge to a stable periodic orbit.Figure 13
Time-series of a trajectory when applying the aPBC withN=62 and (a) k=0 and ko=10; (b) k=0.25 and ko=10; (c) k=0.25 and ko=1000; (d) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)


### 6.1.4. Stabilizing TUPO2 through aPBC: An Odd-Number Limitation Example

One of the main purposes when developing PBC for discrete-time systems is the odd-number limitation subjected by the DFC [14, 21]. Herein, we show a case where the continuous-time DFC with a constant control gain is subjected to the odd-number limitation and aPBC is not. TUPO2 has an odd number of real Floquet multipliers larger than +1 (see Table 1) and is used for a numerical example.Applying the DFC.The DFC with control signal (37)ut,xt=Kxt-T-xt,K∈R2×2, is applied to system (32) with the goal of stabilizing TUPO2. The procedure of finding K for DFC follows [29].A scalar control gaink is used with K defined in (34) resulting in the root locus chart of Figure 14. This chart contains the free system orbit Floquet multipliers μ marked with ×: in black solid line, the fourth largest magnitude (the closed-loop system is infinite-dimensional when applying DFC; see [29] for details) closed-loop Floquet multipliers of the orbit for positive k and, in dark gray, the largest closed-loop Floquet Multiplier of the orbit for negative k. Arrows indicate increasing |k|. It is possible to conclude that for any k there is a Floquet multiplier with magnitude larger than 1 and the orbit is not stabilizable. In addition, all the Floquet multipliers of the orbit tend to +1 when increasing DFC positive scalar gain k. The same chart was done for the other three possible scalars k with equivalent (omitted) results.Figure 14
Root locus chart for the Floquet multipliers of the TUPO2 using the DFC. In black, the four largest magnitude Floquet multipliers for positive control gaink. In dark gray, the largest magnitude Floquet multiplier for negative control gain k. The Floquet multipliers of the free system orbit are represented by (×) and the unit circle by dashed line in black.Following [29], a matrix gain was tried with values obtained by an optimization process using the MATLAB routinesfminsearch andfmincon [35–38]. This optimization aims at finding a minimum for the largest magnitude Floquet multiplier through changing matrix K. A scan process varying all the values of the matrix gain was also tried. The final trial used a matrix gain composed of a rotational matrix (formed by sines and cosines) multiplied by a real scalar gain equivalent to the gain used in [18–20], varying both the rotation angle and the gain. After performing the three methodologies, a matrix control gain that stabilizes TUPO2 was not found.Applying aPBC and Stabilizing TUPO2.Figure 15 shows the magnitude of the Floquet multipliers as a function of the control gain k for the TUPO2 by applying the aPBC. It was noticed that the aPBC stabilizes orbits with an odd number of real Floquet multipliers larger than +1. One detail is relevant here; the stabilization is achieved only for negative values of the gain k and positive values of k make the largest magnitude Floquet multiplier increase.Figure 15
Magnitude of the Floquet multipliers for different values ofk for the TUPO2 controlled by the aPBC.TUPO2 stabilized by the aPBC withk=-0.125 and initial condition x(0)=0.11.8T is shown in state space in Figure 16(a). Time-series of the control signal is presented in Figure 16(b) with u2(t)→0 as t→+∞ indicating stabilization.Figure 16
TUPO2 of forced vdP oscillator stabilized using aPBC withk=-0.125 and x0=[0.11.8]T. (a) State space trajectory and (b) time-series of the control effort.
(a)
(b)


## 6.1.1. Tuning the aPBC Parametersko and N

In order to evaluate the future state estimation, whose characteristics are related to the parametersko and N, we first set k=0. The system is integrated for t∈[0,T] and the points m^j(t+jT/N), j=0,…,N, are collected for different values of ko and N. One then assesses the convergence by comparison with the solution obtained with N=152 and ko=50, by computing (35)εN,ko=m^j,i152,50·-m^j,iN,ko·L10,Tm^j,i152,50·L10,T.The superscript on m^j,i(N,ko)(t) indicates the value of N and ko used in the estimation and the subscript i refers to the state variable of the n-dimensional vector m^j(t).Figure5(a) shows ε(N,50) for different values of N using initial condition x(0)=-0.750.75T,  T=2π/0.45, and state variable i=2. Note that the collocation parameters si and ωj are previously computed for each N tested. The computation of ωj for large N requires intense off-line computational burden and reduces the number of values N tested. These values are computed once for all and are independent of f(t,x(t)). The values of N<72 lead to integration instability and are not computed.Figure 5
Relative error of the free system response estimated future value using different (a)N and (b) ko for the forced vdP oscillator.
(a)
(b)Repeating the same process, we evaluateε(152,ko) for different values of ko; see Figure 5(b). Increasing N and ko results in larger computation effort and trade-off between estimation quality and computation effort should be considered.Due to the error level shown in Figure5 for N=102 and ko=10, we adopt these values in the sequel. Figure 6 shows x2(t) and m^j,2(t-T) (beware of the time shift here). In this case, it is expected that x(t)=φ(t,t-T,x(t-T),0)≈z(t,t-T,x(t-T))≈m^N(102,10)(t-T) and the proximity of both time-series in Figure 6 allows the use of mj(t) instead of m^j(102,10)(t). Once the computational cost of the aPBC comes from processing the set of nN ODEs used for prediction with reasonable error tolerance, m^j(102,10)(t) is also a reference for a good trade-off between future state estimation and computation cost.Figure 6
Current time and shifted future time trajectories of the forced vdP oscillator fork=0, ko=10, and N=102.


## 6.1.2. Target Unstable Periodic Orbits

Two UPOs are chosen as target UPO (TUPO) to be stabilized by the aPBC. TUPO1 was used in [29] for the application of an optimal DFC and in [28] for comparison between aPBC and optimal DFC. Herein, it is used to exemplify the effect of low-precision future state estimation through a small number of collocation points N. TUPO2 is used to exemplify aPBC stabilizing an orbit subjected to the odd-number limitation when applying DFC.Figure7 shows the state space of the target periodic orbits, in continuous-time (line) and Poincaré map (∗), embedded in the chaotic attractor. One point of the orbits with adequate precision for open-loop simulations and their Floquet multipliers are presented in Table 1. Period-1 (where period-1 means one cycle of the forced term or a one-point Poincaré map for the nonautonomous system) orbits of the forced van der Pol oscillator with ω=0.45 leads to T=2π/0.45.Table 1
Floquet multipliers and one Poincaré map point fora=0.988 of the two target periodic orbits of the vdP oscillator.
TUPO1
TUPO2
x
∗
(
t
0
)
[
0.1588445
-
0.1106056
]
T
[
0.0771541
1.6893651
]
T
μ
[−1.872 −0.048]
[
4.627
0.147
]Figure 7
Target UPOs (TUPO1 and TUPO2) embedded in the chaotic attractor (CA) of the forced vdP oscillator. The points for the target UPOs and the chaotic attractor in the Poincaré map (PM) are represented by asterisk and gray bullet, respectively.
(a)
(b)


## 6.1.3. Stabilizing TUPO1 through aPBC for Low-Precision Future State Estimation

In Section6.1.1, the values of N=102 and ko=10 that result in a good estimation of future states according to value ε(N,ko) have been obtained. Results for these N and ko and a comparison with DFC are shown in [28]. In the present section, the aPBC is applied to stabilize TUPO1 using N<102 exemplifying the robustness of the method to the use of low-precision future state estimates, aiming at possible experimental applications.Before applying aPBC to TUPO1, a control gaink that leads to closed-loop stabilization shall be found.Tuning ofk.We use (26), with x∗(0)=0.1588445-0.1106056T, to obtain a value k that stabilizes the orbit of (7a) and (7b). Note that, once the estimation with N=102 and ko=10 is precise, a stabilizing k for (7a) and (7b) yields stability for a slightly perturbed cycle of (32) with the control signal (33) [28].Figure8 shows the magnitude of the two Floquet multipliers |μ| for different values of k. Stability is achieved when |μi|≤1,i=1,2. We choose k=0.25, which results in |μ|max⁡≈0.01213. All the Floquet multipliers in this case are real numbers.Figure 8
Magnitude of the two Floquet multipliers for different values ofk for the TUPO1 controlled by the aPBC.Applying the aPBC for a LowN.Figure 5(a) shows an increase of ε(N,50) when decreasing N for k=0. This value refers to the free system and now we are then interested in characteristics of the closed-loop system, particularly the value of steady state error and control effort as a function of N.Steady state error indexd is the distance (Euclidean norm) between predicted state m^N(t) and current state x(t) for an entire cycle in the close vicinity of target UPO. Steady state control effort index υ is the measure of total external signal (norm-1) necessary to keep the trajectory in close vicinity around one cycle of the target periodic orbit.(36)d=limt→+∞⁡m^Nτ-xτL2t,t+T.υ=limt→+∞⁡uτL1t,t+T.Figure9(a) shows d when stabilizing the TUPO1 by applying the aPBC for different values of N with k=0.25 and ko=10. Figure 9(b) shows the same for υ.Figure 9
Indexes (a)d and (b) υ when applying aPBC to TUPO1 for different N with k=0.25 and ko=10. Bottom dashed line are the values of (a) d and (b) υ computed with N=152 and ko=50 used as reference.
(a)
(b)It can be observed in Figure9 that even for a reduction in the number of collocation points, stability of (an orbit close to) TUPO1 is achieved for N≥67. TUPO1 stabilized with the aPBC for N=82 and N=67 and x(0)=-1.51.0T are shown in Figures 10 and 11, respectively. Figures 10(e) and 11(e) show the initial value for the predicted states mj(0), j=1,2,…,N, and φ(T,0,x(0),0). It is observed that mj(0) is relatively close to φ(T,0,x(0),0).Figure 10
aPBC applied withN=82, k=0.25, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c) and (d)) time-series of the actual and predicted state variables; (e) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)Figure 11
aPBC applied withN=67, k=0.25, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c) and (d)) time-series of the actual and predicted state variables; (e) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)The high values observed in Figure9 for N=3 are explained in Figure 12. It exemplifies a situation wherein a stable orbit different from TUPO1 is achieved due to the reduced N. High amplitude periodic control signal in Figure 12(b) indicates that this stabilized orbit is not an orbit of the original free system. Figure 12(e) shows that N=3 leads to mj(0) almost unrelated to φ(T,0,x(0),0) and consequently u2(t)=k(m^N,2(t)-x2(t)) is a high amplitude control signal. According to Figure 12(a), a trajectory with initial condition x(0) on TUPO1 (Poincaré point) diverges from it showing that it is not stable in this case.Figure 12
aPBC applied withN=3, k=0.25, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c) and (d)) time-series of the actual and predicted state variables; (e) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)It is noticed from the time-series of trajectories in Figures10(c), 10(d), 11(c), and 11(d) that in certain torsions of the trajectory of the system controlled by the aPBC there is an oscillation on m^N(t). This oscillation is generated because small N leads to instability on the predicted state and it reflects in time-series of control signals in Figures 10(b) and 11(b) and on the values of υ and d in Figure 9.Another consequence on reducingN is the loss of stability of m^N(t), exemplified in Figure 13 and occurring even after changing the error tolerances of the explicit R-K integrator (the simulations present in this work were performed using the MATLAB® integration routine ode45 and, in this specific case, the routines ode113 and ode15s were also tried, all of them resulting in interruption of the numerical integration). One of the reasons for this loss of stability is that N is the parameter that rules the implicit R-K integration error and Figure 13(d) shows that mN(0) slightly diverges from φ(T,0,x(0),0) at time T and, contrasting to example of Figure 12, the trajectories of the new system do not converge to a stable periodic orbit.Figure 13
Time-series of a trajectory when applying the aPBC withN=62 and (a) k=0 and ko=10; (b) k=0.25 and ko=10; (c) k=0.25 and ko=1000; (d) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)


## 6.1.4. Stabilizing TUPO2 through aPBC: An Odd-Number Limitation Example

One of the main purposes when developing PBC for discrete-time systems is the odd-number limitation subjected by the DFC [14, 21]. Herein, we show a case where the continuous-time DFC with a constant control gain is subjected to the odd-number limitation and aPBC is not. TUPO2 has an odd number of real Floquet multipliers larger than +1 (see Table 1) and is used for a numerical example.Applying the DFC.The DFC with control signal (37)ut,xt=Kxt-T-xt,K∈R2×2, is applied to system (32) with the goal of stabilizing TUPO2. The procedure of finding K for DFC follows [29].A scalar control gaink is used with K defined in (34) resulting in the root locus chart of Figure 14. This chart contains the free system orbit Floquet multipliers μ marked with ×: in black solid line, the fourth largest magnitude (the closed-loop system is infinite-dimensional when applying DFC; see [29] for details) closed-loop Floquet multipliers of the orbit for positive k and, in dark gray, the largest closed-loop Floquet Multiplier of the orbit for negative k. Arrows indicate increasing |k|. It is possible to conclude that for any k there is a Floquet multiplier with magnitude larger than 1 and the orbit is not stabilizable. In addition, all the Floquet multipliers of the orbit tend to +1 when increasing DFC positive scalar gain k. The same chart was done for the other three possible scalars k with equivalent (omitted) results.Figure 14
Root locus chart for the Floquet multipliers of the TUPO2 using the DFC. In black, the four largest magnitude Floquet multipliers for positive control gaink. In dark gray, the largest magnitude Floquet multiplier for negative control gain k. The Floquet multipliers of the free system orbit are represented by (×) and the unit circle by dashed line in black.Following [29], a matrix gain was tried with values obtained by an optimization process using the MATLAB routinesfminsearch andfmincon [35–38]. This optimization aims at finding a minimum for the largest magnitude Floquet multiplier through changing matrix K. A scan process varying all the values of the matrix gain was also tried. The final trial used a matrix gain composed of a rotational matrix (formed by sines and cosines) multiplied by a real scalar gain equivalent to the gain used in [18–20], varying both the rotation angle and the gain. After performing the three methodologies, a matrix control gain that stabilizes TUPO2 was not found.Applying aPBC and Stabilizing TUPO2.Figure 15 shows the magnitude of the Floquet multipliers as a function of the control gain k for the TUPO2 by applying the aPBC. It was noticed that the aPBC stabilizes orbits with an odd number of real Floquet multipliers larger than +1. One detail is relevant here; the stabilization is achieved only for negative values of the gain k and positive values of k make the largest magnitude Floquet multiplier increase.Figure 15
Magnitude of the Floquet multipliers for different values ofk for the TUPO2 controlled by the aPBC.TUPO2 stabilized by the aPBC withk=-0.125 and initial condition x(0)=0.11.8T is shown in state space in Figure 16(a). Time-series of the control signal is presented in Figure 16(b) with u2(t)→0 as t→+∞ indicating stabilization.Figure 16
TUPO2 of forced vdP oscillator stabilized using aPBC withk=-0.125 and x0=[0.11.8]T. (a) State space trajectory and (b) time-series of the control effort.
(a)
(b)


## 6.2. Applying aPBC to the Rössler System

In this section, we provide numerical results on approximated Prediction-Based Control (aPBC) for the Rössler system using control scheme (21) applied to system (30) according to (32). The goals are (i) to show that the method is applicable to autonomous systems without any modification and (ii) to compare the effect of different numbers of collocation points N on stabilization.aPBC is applied following the same procedure shown in Section6.1 with control signal u(t,x(t)) defined in (33). Herein, we use the constant control gain(38)K=00000000k,k∈R.This choice of matrix K makes only the third state variable influence on feedback and simplifies its tuning.


### 6.2.1. Tuning the aPBC Parametersko and N

The procedure shown in Section6.1.1 is repeated here to characterize the performance of the estimator for different values of the parameters ko and N for k=0. We used initial condition x(0)=0-60.0375T to compute ε(N,ko) in (35). The results are shown in Figure 17 for different values of (a) N with ko=50 (ε(N,50)) and (b) ko with N=152 (ε(152,ko)).Figure 17
Relative error of the free system response estimated future value using different (a)N and (b) ko for the Rössler system.
(a)
(b)In this caseN<32 values lead to numerical integration instability, except when N=3, and are not computed. Herein we notice that less collocation points are necessary to obtain estimation quality equivalent to the one obtained for the vdP oscillator (see Figure 5). It can be seen that for the vdP there is ε(102,50)≈4.4×10-4, while for the Rössler system there is ε(52,50)≈4.4×10-4. Values ε(152,ko) for large ko of the chosen trajectory of the Rössler system also stabilize to a value much smaller than for the vdP oscillator.Figure18 shows x2(t) and m^N,2(t-T) (beware of the time shift here) for the Rössler system, which can be compared with Figure 6. Notice that T=5.8439698764 (the period of the TUPO4) for the Rössler system while T=2π/0.45 for vdP oscillator, but the reduction on ε(152,ko) is much more related to the behavior of the trajectories than the period itself, because collocation points si are computed for a normalized time between 0 and 1.Figure 18
Current time and shifted future time trajectories of the Rössler system fork=0, ko=10, and N=102.


### 6.2.2. Target Unstable Periodic Orbit

Asenjo et al. [34] present one UPO of the chaotic attractor for β1=β2=0.2 and β3=4.5. We apply the Newton-Raphson method with Poincaré section (31) to refine the position and period of the orbit and call it TUPO4 (Figure 19). One point of TUPO4 and its period T with adequate precision for open-loop simulations and Floquet multipliers μ are presented in Table 2.Table 2
Poincaré map pointx(t0), period T, and Floquet multipliers μ of a periodic orbit (TUPO4) of the Rössler system.
TUPO4
x
(
t
0
)
0
5.0206004746
1.6806432235
T
T
5.8439698764
μ
-
1.918
1
-
1.866
×
1
0
-
11Figure 19
TUPO4 of Rössler system in continuous-time (line) and Poincaré map point (∗) for β1=β2=0.2 and β3=4.5.


### 6.2.3. Stabilizing TUPO4 through aPBC

This section is dedicated to the stabilization of TUPO4 using aPBC for different valuesN. From Figure 17, values N=102 and ko=10 are indicated for an accurate estimation of future states and aPBC will be initially applied for these values. Results for low N are also presented, evidencing the robustness of the controller for low-precision future states estimation in autonomous systems.The control gaink of (38) shall be tuned for TUPO4 prior to the application.Tuning ofk.Figure 20 shows the magnitude of the three Floquet multipliers |μ| for different values of k according to the procedure presented in Section 6.1.3. The value μ1=1 is characteristic of autonomous system trajectories and stability is achieved when, for the other Floquet multipliers, |μi|≤1,i=2,3. Choosing k=1.85 results in |μ2|≈2.994×10-3 (and |μ3|≈2.167×10-14). All Floquet multipliers in this case are real numbers.Figure 20
Magnitude of the Floquet multipliers for different values ofk for the TUPO4 controlled by the aPBC.Applying the aPBC for Different Values ofN.In this section, we evaluate the characteristics of aPBC for different values N applied to the stabilization of TUPO4. Herein, we use N≤102,  ko=10, and k=1.85.Figure17(a) shows an increase of ε(N,50) observed when decreasing N for k=0. Figure 21 shows d and υ for different values of N with k=1.85 and ko=10. Values of 3<N<32 were not computed because of instability, the same as observed in Section 6.1.3.Figure 21
Indexes (a)d and (b) υ when applying aPBC to TUPO4 for different N with k=1.85 and ko=10. The bottom dashed lines are the values of (a) d and (b) υ computed with N=152 and ko=50 used as reference.
(a)
(b)Herein, we observe thatN≥102 leads to a precise estimation once there is a convergence of υ to the value computed using N=152 and ko=50 (dashed line). Comparing Figure 21 with Figure 9, a great reduction on υ and d is observed for the control method applied to the Rössler system, which is in accordance with what was observed for ε(N,50). In fact, N between 47 and 52 leads to a value of d equivalent to the value obtained for N=102 for the vdP oscillator. Considering υ, this reduction is more significant, with N=37 leading to a value of υ close to the value obtained for N=102 for the vdP. This result is important for practical aspects because it shows that, for certain systems, the augmented state of the closed-loop system controlled with the aPBC can be reduced significantly (in comparison with the initial results shown for the vdP oscillator).Figure21 shows that stability of TUPO4 is achieved for N≥32. The TUPO4 stabilized by the aPBC for N=102, N=52, and N=32 is shown in Figures 22, 23, and 24, respectively, using x(0)=0-60.0375T. Deviation of m^N,2(t) from x2(t) is observed in Figure 24(d), which reflects the numerical integration instabilities for N<32. We verify that stabilization is not achieved for N=3, which leads to the high values of υ and d observed in Figure 21.Figure 22
aPBC applied to the Rössler system withN=102, k=1.85, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c), (d), and (e)) time-series of the actual and predicted state variables; (f) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)
(f)Figure 23
aPBC applied to the Rössler system withN=52, k=1.85, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c), (d), and (e)) time-series of the actual and predicted state variables; (f) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)
(f)Figure 24
aPBC applied to the Rössler system withN=32, k=1.85, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c), (d), and (e)) time-series of the actual and predicted state variables; (f) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)
(f)


## 6.2.1. Tuning the aPBC Parametersko and N

The procedure shown in Section6.1.1 is repeated here to characterize the performance of the estimator for different values of the parameters ko and N for k=0. We used initial condition x(0)=0-60.0375T to compute ε(N,ko) in (35). The results are shown in Figure 17 for different values of (a) N with ko=50 (ε(N,50)) and (b) ko with N=152 (ε(152,ko)).Figure 17
Relative error of the free system response estimated future value using different (a)N and (b) ko for the Rössler system.
(a)
(b)In this caseN<32 values lead to numerical integration instability, except when N=3, and are not computed. Herein we notice that less collocation points are necessary to obtain estimation quality equivalent to the one obtained for the vdP oscillator (see Figure 5). It can be seen that for the vdP there is ε(102,50)≈4.4×10-4, while for the Rössler system there is ε(52,50)≈4.4×10-4. Values ε(152,ko) for large ko of the chosen trajectory of the Rössler system also stabilize to a value much smaller than for the vdP oscillator.Figure18 shows x2(t) and m^N,2(t-T) (beware of the time shift here) for the Rössler system, which can be compared with Figure 6. Notice that T=5.8439698764 (the period of the TUPO4) for the Rössler system while T=2π/0.45 for vdP oscillator, but the reduction on ε(152,ko) is much more related to the behavior of the trajectories than the period itself, because collocation points si are computed for a normalized time between 0 and 1.Figure 18
Current time and shifted future time trajectories of the Rössler system fork=0, ko=10, and N=102.


## 6.2.2. Target Unstable Periodic Orbit

Asenjo et al. [34] present one UPO of the chaotic attractor for β1=β2=0.2 and β3=4.5. We apply the Newton-Raphson method with Poincaré section (31) to refine the position and period of the orbit and call it TUPO4 (Figure 19). One point of TUPO4 and its period T with adequate precision for open-loop simulations and Floquet multipliers μ are presented in Table 2.Table 2
Poincaré map pointx(t0), period T, and Floquet multipliers μ of a periodic orbit (TUPO4) of the Rössler system.
TUPO4
x
(
t
0
)
0
5.0206004746
1.6806432235
T
T
5.8439698764
μ
-
1.918
1
-
1.866
×
1
0
-
11Figure 19
TUPO4 of Rössler system in continuous-time (line) and Poincaré map point (∗) for β1=β2=0.2 and β3=4.5.


## 6.2.3. Stabilizing TUPO4 through aPBC

This section is dedicated to the stabilization of TUPO4 using aPBC for different valuesN. From Figure 17, values N=102 and ko=10 are indicated for an accurate estimation of future states and aPBC will be initially applied for these values. Results for low N are also presented, evidencing the robustness of the controller for low-precision future states estimation in autonomous systems.The control gaink of (38) shall be tuned for TUPO4 prior to the application.Tuning ofk.Figure 20 shows the magnitude of the three Floquet multipliers |μ| for different values of k according to the procedure presented in Section 6.1.3. The value μ1=1 is characteristic of autonomous system trajectories and stability is achieved when, for the other Floquet multipliers, |μi|≤1,i=2,3. Choosing k=1.85 results in |μ2|≈2.994×10-3 (and |μ3|≈2.167×10-14). All Floquet multipliers in this case are real numbers.Figure 20
Magnitude of the Floquet multipliers for different values ofk for the TUPO4 controlled by the aPBC.Applying the aPBC for Different Values ofN.In this section, we evaluate the characteristics of aPBC for different values N applied to the stabilization of TUPO4. Herein, we use N≤102,  ko=10, and k=1.85.Figure17(a) shows an increase of ε(N,50) observed when decreasing N for k=0. Figure 21 shows d and υ for different values of N with k=1.85 and ko=10. Values of 3<N<32 were not computed because of instability, the same as observed in Section 6.1.3.Figure 21
Indexes (a)d and (b) υ when applying aPBC to TUPO4 for different N with k=1.85 and ko=10. The bottom dashed lines are the values of (a) d and (b) υ computed with N=152 and ko=50 used as reference.
(a)
(b)Herein, we observe thatN≥102 leads to a precise estimation once there is a convergence of υ to the value computed using N=152 and ko=50 (dashed line). Comparing Figure 21 with Figure 9, a great reduction on υ and d is observed for the control method applied to the Rössler system, which is in accordance with what was observed for ε(N,50). In fact, N between 47 and 52 leads to a value of d equivalent to the value obtained for N=102 for the vdP oscillator. Considering υ, this reduction is more significant, with N=37 leading to a value of υ close to the value obtained for N=102 for the vdP. This result is important for practical aspects because it shows that, for certain systems, the augmented state of the closed-loop system controlled with the aPBC can be reduced significantly (in comparison with the initial results shown for the vdP oscillator).Figure21 shows that stability of TUPO4 is achieved for N≥32. The TUPO4 stabilized by the aPBC for N=102, N=52, and N=32 is shown in Figures 22, 23, and 24, respectively, using x(0)=0-60.0375T. Deviation of m^N,2(t) from x2(t) is observed in Figure 24(d), which reflects the numerical integration instabilities for N<32. We verify that stabilization is not achieved for N=3, which leads to the high values of υ and d observed in Figure 21.Figure 22
aPBC applied to the Rössler system withN=102, k=1.85, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c), (d), and (e)) time-series of the actual and predicted state variables; (f) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)
(f)Figure 23
aPBC applied to the Rössler system withN=52, k=1.85, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c), (d), and (e)) time-series of the actual and predicted state variables; (f) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)
(f)Figure 24
aPBC applied to the Rössler system withN=32, k=1.85, and ko=10. (a) Trajectory in state space; (b) time-series of the control signal; ((c), (d), and (e)) time-series of the actual and predicted state variables; (f) time-series of φ(T,0,x(0),0) and the initial value of the predicted states mj(0).
(a)
(b)
(c)
(d)
(e)
(f)


## 7. Conclusions

Approximate Prediction-Based Control (aPBC), which is the continuous-time version of Prediction-Based Control (PBC), was reviewed in this work and some of its characteristics were analyzed through numerical simulations. It was shown that continuous-time application of PBC requires future state values and aPBC is based on a methodology for estimating them in real-time.The main drawback in aPBC is that estimation of future states requires extending the order of the closed-loop system by the number of implicit Runge-Kutta solution points, herein the orthogonal collocation points, times the number of state variables. Applications require real-time explicit integration of these extended equations and reducing them implies reducing the computational cost while reducing estimation precision. Both exemplified cases have shown that there is a lower bound on the number of collocation points to achieve stability, to both target UPO and numerical integration. While comparing examples, it was noticed that the number of collocation points depends upon the trajectory torsions in state space. Future studies should be carried out in order to reduce the number of implicit Runge-Kutta solution points through different orthogonal polynomials or different implementation methods.Another interesting example is an UPO of a continuous-time system that could not be stabilized by DFC with constant control gain because of the odd-number limitation. The same orbit is stabilized by aPBC with a negative scalar constant gain without any modification on the method. Then, the main justification for the application of PBC in comparison to DFC in discrete-time is also valid for continuous-time.The method proposed for estimation of open-loop future states of nonlinear systems applied to aPBC may be generalized for other applications. As example, model-based predictive control [39] applied to nonlinear systems requires open-loop future state values to optimize control gain at each time step and this routine may be fastened by using the proposed scheme.



---

*Source: 3298286-2018-03-01.xml*