Adv Syst Sci Appl 2020; 01:27–49
Published online at http://ijassa.ipu.ru/index.php/ijassa/article/view/775

Average Predictive Control for Nonlinear Discrete Dynamical
Systems

Dmitriy Dmitrishin1, Emil Iacob2∗, Alex Stokolos2

1Odessa National Polytechnic University, Odessa, Ukraine
2Georgia Southern University, Statesboro GA 30460, USA

Abstract: We explore the problem of stabilization of unstable periodic orbits in discrete nonlinear
dynamical systems. This work proposes the generalization of predictive control method for
resolving the stabilization problem. Our method embodies the development of control method
proposed by B.T. Polyak. The control we propose uses a linear (convex) combination of iterated
functions. With the proposed method auxiliary, the problem of robust cycle stabilization for
various cases of its multipliers localization is solved. An algorithm for finding a given length
cycle when its multipliers are known is described as a particular case of our method application.
Also, we present numerical simulation results for some well-known mappings and the possibility
of further generalization of this method.

Keywords: Non-linear discrete systems, stabilization, predictive control.

1. INTRODUCTION

Nonlinear dynamical systems are often characterized by extremely unstable movements in
the phase space, defined as chaotic movements [5]. In practice, it is generally desirable to
suppress or prevent such chaotic behavior due to its adverse effect on the physical systems
normal operation. Due to its theoretical significance and engineering applicability, much
attention has been paid to the problem of chaos controlling in various fields and numerous
studies [2, 3, 15, 22]. By chaos control we mean small external influences on the system or a
small change in the system structure in order to transform the system chaotic behavior into a
regular (or chaotic, still specific with other properties) one [13].

It is assumed that the dynamical system includes a chaotic attractor, which contains
a countable set of unstable in Lyapunov sense cycles with different periods. If, by using
some control effect, a certain cycle is stabilized locally, the system path will remain in its
neighborhood, i.e. regular movements will be observed in the system. Hence, one of the ways
for chaos controlling refers to the local stabilization of certain orbits from a chaotic attractor.

The problem of stabilizing cycles is closely related to the problem of finding periodic
points. The various control schemes [1, 4, 9] that were proposed for solving these problems

∗Corresponding author: ieiacob@georgiasouthern.edu



28 D. DMITRISHIN, E. IACOB, A. STOKOLOS

can be divided into two large groups: direct and indirect methods. The indirect methods either
use the initial mapping T iterations or imply the construction a system of which order is
T times greater than the initial system order (T being the desired cycle length). Then one
of the methods of finding a fixed point is applied. The most common among fixed point
finding techniques is the Newton-Raphson relaxation method and its further modifications
[12, 19, 22, 31]. The next step is to select periodic points from the entire set of fixed points.
In direct methods, all the points in the cycle are found concurrently, i.e. the whole cycle
is stabilized. In this case, the initial system is closed by control, based on the feedback
principle [7, 20, 25, 30]. Among such control schemes, the most simple in terms of physical
implementation are the linear ones. However, they have significant limitations, as they can
only be applied to a narrow domain of the space of parameters that are part of the initial
nonlinear system [6, 28, 32]. To overcome the restricted stabilization condition Ushio and
Yamamoto [29] introduced a prediction-based feedback control method for discrete chaotic
systems with accurate mathematical model. B.T. Polyak [23] (see also [24]) proposed a direct
predictive control method that works well for one-dimensional as well as multidimensional
maps. L. Shalby [27] used predictive feedback control method for stabilization of continuous
time systems.

Another possible control schemes classification into two groups: methods using the Jacobi
matrix and methods not based on this matrix. Naturally, it is assumed that the Jacobi matrix
at the cycle points is not properly known, otherwise it would be possible to use the whole
powerful apparatus of the linear control theory applied to systems linearized in the cycle
neighborhood. The Jacobi matrix is an indispensable attribute of Newton-Raphson-type
methods. This matrix is also used in one of the modifications derived from Polyak predictive
control method.

One of the main disadvantages of Polyak scheme refers to the need for knowing the Jacobi
matrix for the cycle, or at least the need for sufficiently good estimates of the cycle multipliers.
The research exposed herein is purposed to improve the Polyak method by replacing it with
mixed predicted values. A Jacobi matrix representation for the T cycle at a controlled system
is found through the Jacobi matrix of the same cycle in the initial system. Therefore, the
correspondence between the cycle multipliers of the open loop system and those of closed
loop system is established.

Below, it is assumed that the initial system cycle multipliers are not exactly known, we
only know the range of their localization. Then the solution of cycles robust stabilization
problem for various localization of multipliers is given, and taking into account these general
provisions the Polyak method is considered. It is worth noting that in the general case
of complex multipliers, we must know precisely enough their localization regions. At the
end, the applications of proposed predictive control scheme to stabilize the cycles of some
common systems in Physics literature are considered.

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



AVERAGE PREDICTIVE CONTROL 29

2. PROBLEM STATEMENT

Considered is a nonlinear discrete system

xn+1 = f (xn) , xn ∈ Rm, n = 1, 2, . . . , (2.1)

where f(x) is a differentiable vector function of corresponding dimension. The
differentiability is redundant and can be relaxed (f(x) needs to be differentiable only at
the cycle points, c.f. [23]). It is assumed that the system (2.1) has an invariant convex set
A, i.e. if ξ ∈ A, then f(ξ) ∈ A. We emphasize that we do not assume that the set A is a
minimally convex set. It is also assumed that this system has one or several unstable T
cycles {η1, . . . , ηT}, where all vectors are different and belong to the invariant set A, i.e.
ηj+1 = f (ηj), j = 1, . . . , T − 1, η1 = f (ηT ). The considered unstable cycles multipliers are

determined as eigenvalues of the product of Jacobi matrices
T∏
j=1

f ′ (ηT−j+1) of dimensions

m×m at the cycle points. The matrix
T∏
j=1

f ′ (ηT−j+1) is called the Jacobi matrix of the

cycle {η1, . . . , ηT} . Typically, a priori the cycles of the system (2.1) are not known.

Consequently, the spectrum {µ1, . . . , µm} of the matrix
T∏
j=1

f ′ (ηT−j+1) is unknown as well.

The spectrum elements are called cycle multipliers. Below, we assume that some estimates
on the localization set M for the cycle multipliers are known.

Let us consider the control system

xn+1 = f(xn) + un, (2.2)

where un = (θ1 − 1)f(xn) +
N∑
j=2

θjf
((j−1)T+1)(xn), f (1)(x) = f(x), f (k)(x) =

f
(
f (k−1)(x)

)
, k = 2, . . . , (N − 1)T + 1. The numbers θ1, . . . , θN are real. It can be

easily verified that at
N∑
j=1

θj = 1 the system (2.2) also includes the cycle {η1, . . . , ηT}. We

aim to choose such parameter N and coefficients {θ1, . . . θN} so that the system (2.2) cycle
{η1, . . . , ηT} would be locally asymptotically stable. Naturally, when constructing these
coefficients, there will be used information on set M of multipliers localization. It is also
desirable [16, 17] to fulfill an additional condition: the system (2.1) invariant convex set A
must be also invariant for the system (2.2). This requirement will be fulfilled, for example, if
0 ≤ θj ≤ 1, j = 1, . . . , N . )

Polyak method [23] utilizes the case θ1 = 1, θ2 = · · · = θN−2 = 0, θN−1 = −θN = ε.
Regarding the set M , it was assumed that M = D ∪ {µ∗} where D = {z : |z| < 1} is the
central unit circle on the complex plane, and µ∗ is a known real number. In the case m = 1

the required coefficient formula has the form ε =
1∓ (|µ∗| /ρ)−

1
T

(µ∗)N−2 (µ∗ − 1)
where 0 < ρ < 1. In this

article the control problem is solved for a wider class of multipliers localization set M .

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



30 D. DMITRISHIN, E. IACOB, A. STOKOLOS

3. CONSTRUCTING THE JACOBI MATRIX FOR A CONTROLLED SYSTEM

Investigating stability of T cycles of the system (2.2) consists in constructing of Jacobi matrix
T∏
j=1

F ′ (ηT−j+1) of that cycle and studying the eigenvalues of this matrix. To derive the Jacobi

matrix, we use the ideas from [23].
Let Jj = f ′ (ηj), j = 1, . . . , T , then we write the Jacobi matrix of the system (2.1) cycle

{η1, . . . , ηT} as J = JT · . . . · J1. We introduce the following auxiliary matrices:

A1 = I, A2 = J1, A3 = J2 · J1, . . . , AT−1 = JT−1 · . . . · J1

(I − unity matrix of order m×m);

B1 = JT · . . . · J1 = J, B2 = JT · . . . · J2, . . . , BT = JT ;

then BkAk = J , k = 1, . . . , T , AkBk = (Jk−1 · . . . · J1) · (JT · . . . · Jk) and, consequently,
(AkBk)

s = AkJ
s−1Bk, s = 1, 2, . . ..

By chain rule:(
f (s)(x)

)′∣∣∣
x=ηi

=
(
f (s−1)(x)

)′∣∣∣
x=ηi+1

· (f(x))′
∣∣
x=ηi

=
(
f (s−1)(x)

)′∣∣∣
x=ηi+1

· Ji,

we get (
f ((j−1)T )(x)

)′∣∣∣
x=ηi

= AiJ
j−2Bi, j = 2, . . . , N

and therefore (
f ((j−1)T+1)(x)

)′∣∣∣
x=ηi

= JiAiJ
j−2Bi, j = 2, . . . , N.

Next we find that

F ′ (ηi) =
N∑
j=1

θj
(
f ((j−1)T+1)(x)

)′∣∣∣
x=ηi

= θ1Ji +
N∑
j=2

θjJiAiJ
j−2Bi.

For the Jacobi matrix of the system (2.2) cycle {η1, . . . , ηT} we can write:

F ′ (ηT ) · . . . · F ′ (η1) = JT

(
θ1I + AT

(
N∑
j=2

θjJ
j−2

)
BT

)
·

JT−1

(
θ1I + AT−1

(
N∑
j=2

θjJ
j−2

)
BT−1

)
· . . . · J1

(
θ1I + A1

(
N∑
j=2

θjJ
j−2

)
B1

)

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



AVERAGE PREDICTIVE CONTROL 31

Taking into account that

Jk

(
θ1I + Ak

(
N∑
j=2

θjJ
j−2

)
Bk

)
=

JkAk

(
θ1I +

(
N∑
j=2

θjJ
j−2

)
BkAk

)
A−1k = JkAk

(
N∑
j=1

θjJ
j−1

)
A−1k

and JkAk = Ak+1 it follows that

F ′ (ηT ) · . . . · F ′ (η1) =

= JTAT

(
N∑
j=1

θjJ
j−1

)
A−1T · JT−1AT−1

(
N∑
j=1

θjJ
j−1

)
A−1T−1 · . . . ·

· J1A1

(
N∑
j=1

θjJ
j−1

)
A−11 = J

(
N∑
j=1

θjJ
j−1

)T

(The reader is gently advised that the superscript T in the formula above and all
subsequent formulas denotes power, not transpose.)

For deriving the Jacobian formula above it was assumed that the matrix J was not
degenerated. This limitation can be easily circumvented using a well-known topological
technique: considering the matrix J + δI instead of the degenerated matrix J and after all
calculations taking the limit as δ → 0. Thus, the following result is obtained.

Lemma 3.1:
The Jacobi matrix of the cycle {η1, . . . , ηT} in the system (2.2) can be represented as

J

(
N∑
j=1

θjJ
j−1

)T

, (3.3)

where J is the Jacobi matrix of the cycle {η1, . . . , ηT} in the system (2.1).

We now consider another control system, instead of system (2.2):

xn+1 = f

(
θ1xn +

N∑
j=2

θjf
((j−1)T ) (xn)

)
. (3.4)

When
N∑
j=1

θj = 1 then the system (3.4) preserves the cycle {η1, . . . , ηT} . In addition,

according to formula (3.3), the Jacobi matrix of the system (3.4) cycle is expressed in the
terms of Jacobi matrix of the system (2.1). The advantage of the control system (3.4) over the
system (2.2) consists of a fewer calculation of the values for function f(x) (more precisely,
the difference is N − 2).

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



32 D. DMITRISHIN, E. IACOB, A. STOKOLOS

4. MAIN RESULT

All results presented in this section are formulated for system (2.2), however they hold
without change for system (3.4).

Theorem 4.1:
Suppose f ∈ C1 and that the system (2.1) has an unstable T cycle with multipliers
{µ1, . . . , µm}. Then this cycle will be a locally asymptotically stable cycle of the system (2.2)
if

µj [r (µj)]
T ∈ D, j = 1, . . . ,m,

where r(µ) =
N∑
j=1

θjµ
j−1.

Proof
According to Lemma 3.1, the characteristic polynomial for a system of linear approximation
in the cycle neighbourhood in the case of system (2.2) can be written as ϕ(λ) =
det
(
λI − J [r(J)]T

)
. By reducing the matrix J to the Jordan form, this characteristic

polynomial can be represented as ϕ(λ) =
m∏
j=1

(
λI − µj [r (µj)]

T
)

, from where the theorem

conclusion follows.

Note that the condition r(1) = 1 is obligatory. If, additionally, θj ∈ [0, 1] for j =

1, . . . , N , then µj [r (µj)]
T ∈ D when µj ∈ D, and hence

∣∣∣µj [r (µj)]
T
∣∣∣ < |µj|1+T . This means

that if some multiplier of the system (2.1) cycle lies in the unit circle, the corresponding
multiplier of the system (2.2) will lie closer to zero. Thus for the closed loop system, the
stabilization quality is improving. Various estimates for multipliers allow us to construct
control systems that stabilize cycles.

4.1. Case M = {µ1, . . . , µm}

If the multipliers are exactly known, we can choose N = m+ 1 and the coefficients

{θ1, . . . , θm+1} from the condition r(µ) =
m+1∑
j=1

θjµ
j−1 = 1

m∏
k=1

(1−µk)

m∏
k=1

(µ− µk). Then from

Theorem 4.1 we get the following conclusion.
Conclusion. Suppose that f ∈ C1 and the system (2.1) has an unstable T cycle with

multipliers {µ1, . . . , µm}, and the coefficients θ1, . . . , θm+1 are found as exposed above. Then
this cycle will be a locally asymptotically stable cycle of system (2.2). Moreover, if the initial
point belongs to the cycle basin of attraction, the convergence to the cycle is superlinear.

The superlinearity of the convergence rate follows from the fact that all multipliers of
system (2.2) {η1, . . . , ηT} cycle turn out to be zero.

Note that the authors are unaware about any other method that allow to stabilize a cycle
by knowing the cycle multipliers only. Unfortunately, in a typical situation the multipliers are
either unknown. The best we can expect is to localize them approximately. What to do in that
case is considered in the next section.

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



AVERAGE PREDICTIVE CONTROL 33

4.2. Case M = {z : Re z ≤ 0} ∪ D

Theorem 4.2:
Suppose f ∈ C1 and that system (2.1) has an unstable T -cycle with multipliers {µ1, . . . , µm}
satisfying the conditions:

|µj − µ̂j| < δj, Re {µj} ≤ 0, j = 1, . . . , n1, |µj| < 1, j = n1 + 1, . . . ,m.

Let the coefficients θj , j = 1, . . . , N , of the system (2.2) be determined from the condition

n1+1∑
j=1

θjµ
j−1 =

1
n1∏
k=1

(1− µ̂k)

n1∏
k=1

(µ− µ̂k) (here N = n1 + 1).

Then, for sufficiently small values δj , j = 1, . . . , n1, the T -cycle will be a locally
asymptotically stable cycle of system (2.2).

Proof
Since Re {µj} < 0, j = 1, . . . , n1, then all coefficients θj > 0, j = 1, . . . , n1. That means

that
∣∣∣µ [r (µ)]T

∣∣∣ < |µ|1+T when |µ| < 1, i.e. the eigenvalues of the Jacobi matrix of
system (2.2) cycle corresponding to multipliers µj , j = n1 + 1, . . . ,m, are smaller than the
multipliers absolute values. Let δj = 0, j = 1, . . . , n1, then the eigenvalues corresponding
to multipliers µj , j = 1, . . . , n1, are equal to zero. When δj , j = 1, . . . , n1, are sufficiently
small in magnitude, these eigenvalues will lie in the central unit circle, as follows from
Rouche theorem. Thus, all eigenvalues will less than 1 in absolute value, which means local
asymptotic stability.

4.3. Case M = D ∪ {µ∗} , |µ∗| > 1 [23]

In [23], the coefficients θ1, . . . , θN were chosen as θ1 = 1, θ2 = . . . = θN−2 = 0, θN−1 =

−θN = ε, where ε =
1∓ (|µ∗| /ρ)−

1
T

(µ∗)N−2 (µ∗ − 1)
, 0 < ρ < 1. Such a choice ensures that the multipliers

belong to the open central unit interval corresponding to µ∗. However, the condition∣∣∣µ [r(µ)]T
∣∣∣ < 1 with |µ| < 1 is not necessarily satisfied. Nevertheless, the value ε can be

made arbitrarily small by choosing the number N large. And then, from Rouche’s theorem,
it follows that with a sufficiently large N the other multipliers will remain within the central
unit circle. This ensures the local asymptotic stability of the system (2.2) cycle.

When it is known that µ∗ < −1, the control scheme can be simplified, namely:

xn+1 = θ1f (xn) + θ2f
(T+1) (xn) , (4.5)

where θ1 =
|µ∗|

1 + |µ∗|
, θ2 =

1

1 + |µ∗|
.

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



34 D. DMITRISHIN, E. IACOB, A. STOKOLOS

4.4. Case M = D ∪ {µ∗, µ∗} , |µ∗| > 1

The case of general localization of multipliers {µ1, . . . , µm} for the system (2.1) cycle was
considered in [23] but only for T = 1. In that case, the coefficients θ1, . . . , θN were no longer
scalars but matrices and, as before, were chosen as θ1 = I , θ2 = . . . = θN−2 = 0, θN−1 =
−θN = ε, where I is identity matrix, 0 is zero matrix, ε = SΛS−1, Λ = diag {ε1, . . . , εm},

εj =
1 + eıϕ (ρ/ |µj|)
(µj)

N−2 (µj − 1)
, if |µj| < 1, and εj = 0, if |µj| > 1, and 0 < ρ < 1, ϕ ∈ {0, π} if

µj as a real number. The matrix S consists of the eigenvectors of the Jacobi matrix J for
equilibrium point. Thus, to apply the stabilization method, it is necessary to know not only
all the multipliers of the equilibrium, but also the Jacobi matrix itself. That is impossible
when the equilibrium is not known.

Now, let us we apply the scheme (2.2). Let µ∗ = ρeıϕ. Then

r(µ) =
(µ− µ∗) (µ− µ∗)
(1− µ∗) (1− µ∗)

=

ρ2

ρ2 − 2ρ cosϕ+ 1
+

−2ρ cosϕ

ρ2 − 2ρ cosϕ+ 1
µ+

1

ρ2 − 2ρ cosϕ+ 1
µ2.

If the complex number µ∗ lies in the left half-plane, then the coefficients of polynomial
r(µ) are positive, so

∣∣∣µ [r(µ)]T
∣∣∣ < |µ|1+T with |µ| < 1. Therefore, each multiplier of system

(2.2) cycle lying in the central unit circle turns out to be in absolute value less then the
multiplier of the corresponding cycle of the system (2.1). Also, the multipliers corresponding
to µ∗ and µ∗ change to zero. If the multipliers µ∗ and µ∗ are not exactly known, but they can
be well estimated, then for the coefficients θ1, θ2, θ3, although different from the calculated
ones, the values of the polynomial with these coefficients at points µ∗ and µ∗ will not exceed
1 in absolute value, as follows from the Rouche theorem. The desired control system is

xn+1 = θ1f (xn) + θ2f
(T+1) (xn) + θ3f

(2T+1) (xn) ,

where θ1 =
ρ2

ρ2 − 2ρ cosϕ+ 1
, θ2 =

−2ρ cosϕ

ρ2 − 2ρ cosϕ+ 1
, θ3 =

1

ρ2 − 2ρ cosϕ+ 1
.

4.5. Case T = 1, M = b−µ∗, µ∗c, M = b−µ∗, 1c

Suppose M = b−µ∗, µ∗c. From Theorem 4.1 it follows that in order to stabilize the
equilibrium, it would be necessary to construct a polynomial µr(µ), so that r(1) = 1 and
|µr(µ)| ≤ 1 for all |µ| < µ∗.

Theorem 4.3:
Let f ∈ C1 and the system (2.1) has unstable equilibrium with multipliers {µ1, . . . , µm} ⊂
[−µ∗, µ∗]. Let the value N be odd and be chosen from the condition csc

π

2N
> µ∗, and the

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



AVERAGE PREDICTIVE CONTROL 35

coefficients θ1, . . . , θN from the condition

µr(µ) = µ

N∑
j=1

θjµ
j−1 = (−1)

N−1
2 TN

(
µ sin

π

2N

)
,

where TN(x) is the first kind Chebyshev polynomial of odd order N . Then this equilibrium
will be a locally asymptotically stable equilibrium of system (2.2) (modulo a finite number of

cases when µj =
cos πk/N

sinπ/2N
for some k = 1, . . . , N − 1).

The proof follows from the properties of the first kind Chebyshev polynomials: |µr(µ)| ≤
1 at

∣∣∣µ sin
π

2N

∣∣∣ ≤ 1, r(1) = 1.Note that µ∗ →∞ (N →∞) with asymptotics
2

π
N .

Now we will consider the case of M = b−µ∗, 1c.

Theorem 4.4:
Let f ∈ C1 and the system (2.1) has unstable equilibrium with multipliers {µ1, . . . , µm} ⊂
b−µ∗, 1c. Let the N value be chosen from the condition cot2 π

4N
> µ∗, and the coefficients

θ1, . . . , θN from the conditions

µr(µ) = µ
N∑
j=1

θjµ
j−1 = TN

(
µ
(

1− cos
π

2N

)
+ cos

π

2N

)
,

where TN(x) is the first kind Chebyshev polynomial of order N . Then this equilibrium will
be locally asymptotically stable equilibrium of the system (2.2) (modulo a finite number of
cases).

Proof
Note that r(0) = TN(cos π/2N) = 0, r(1) = TN(1) = 1. In addition |TN(x)| ≤ 1, at |x| ≤ 1,
whence |µr(µ)| ≤ 1 at

∣∣∣µ(1− cos
π

2N

)
+ cos

π

2N

∣∣∣ ≤ 1. The last inequality is equivalent to

− cot2
π

4N
≤ µ ≤ 1, which proves the theorem.

Note that µ∗ →∞ (N →∞) with asymptotics
16

π2
N2.

4.6. The general case

Using the ideas from Theorem 4.1 cases, we can propose the following T -cycle stabilization
scheme, for which the coefficients θj are not necessarily constants:

a) find the matrix f ′(x),
b) find the vectors f (s)(x), s = 1, . . . , T − 1,
c) find the matrix f ′

(
f (T−1)(x)

)
· . . . · f ′(f(x)) · f ′(x)

d) find the matrix characteristic polynomial
m+1∑
j=1

θj(x)µj−1,

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



36 D. DMITRISHIN, E. IACOB, A. STOKOLOS

e) normalize the characteristic polynomial
1

m+1∑
j=1

θj(x)

m+1∑
j=1

θj(x)µj−1,

f) build the control system
xn+1 = F (xn) ,

where

F (x) =
1

m+1∑
j=1

θj(x)

m+1∑
j=1

θj(x)f ((j−1)T+1)(x)

or

F (x) = f

 1
m+1∑
j=1

θj(x)

(
θ1x+

m+1∑
j=2

θj(x)f ((j−1)T )(x)

)
Let us consider how this scheme looks like in the case of a linear problem. Let f(x) = Ax

whereA is a non-degeneratem×mmatrix. Then η = 0 is a single fixed point, in the absence

of any higher order cycles. We choose θj from the condition
1

det (I − A)
det (µI − A) =

m+1∑
j=1

θjµ
j−1. Then the control system is xn+1 =

1
m+1∑
j=1

θj

m+1∑
j=1

θjA
jxn. By the Hamilton-Cayley

theorem it follows that this system right-hand side is an identical zero.
In the general case, applying this method to stabilizing chaotic motion tending to mixing,

one can expect that after a certain number of iterations the trajectory falls into the basin of
attraction for the stabilized cycle. Then the convergence to the cycle will be superlinear.

Note that if in all the considered cases |θj| is being used instead of θj , then it becomes
possible to stabilize the system (2.1) cycles with multipliers lying inM = D ∪ {µ : Re (µ) ≤
0}. Moreover, the convex invariant set of system (2.1) will remain such for system (2.2). In
addition, the system (2.2) multipliers, corresponding to those multipliers of system (2.1) that
lie in the unit circle, will become closer to zero.

5. EXAMPLES

Let us illustrate the effectiveness of the generalized predictive control method for finding
periodic orbits with several well-known examples of scalar and vector chaotic systems [26].

The scheme applied for the logistic and triangular mappings (examples from sections
5.1–5.2) was a general scheme{

xn+1 = θ(xn)
1+θ(xn)

f (xn) + 1
1+θ(xn)

f (T+1) (xn) ,

θ (xn) = −f ′
(
f (T−1) (xn)

)
· . . . · f ′ (f (xn)) · f ′ (xn) .

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



AVERAGE PREDICTIVE CONTROL 37

It was possible to find a large number of cycles for all considered periods T ; in general,
different initial conditions are producing different cycles. Numerical calculations show that
with sufficiently dense initial values grid, all cycles of a given length can be found. However,
in this case it is necessary to ensure that the point xn remains within the invariant set,
otherwise, as a rule it goes to infinity. If we use |θ(x)| instead of θ(x), the point xn will always
remain in the invariant set. However, in this case we can find cycles only with multipliers from
the set M = D ∪ {µ : Re (µ) ≤ 0}.

In the two-dimensional case, the scheme used was

xn+1 = f

(
θ

1 + θ
xn +

1

1 + θ
f (T ) (xn)

)
. (5.6)

The value θ should be chosen according to the condition

x

(
θ

1 + θ
+

1

1 + θ
x

)T
∈ D

at x = µ∗j , where µ∗j are cycle multipliers (j = 1, 2), and in general, they are unknown. In the
examples below, one of the multipliers never exceeds one in magnitude, while the second one
is negative, greater than one in absolute value.

The Theorem 4.2 guarantees stability conditions even if the parameter θ satisfies

x

(
θ

1 + θ
+

1

1 + θ
x

)T
6= 0.

It is enough have θ in the neighborhood of multiplier. In general Theorem 4.2 does not provide
the estimate on the parameters, however when one multiplier is in the unit disc and the other
is real and negative an elementary trial works quite effectively.

Thus θ > 0, and we only have to check the compliance with the condition for the second
multiplier ∣∣∣∣∣µ∗2

(
θ

1 + θ
+

1

1 + θ
µ∗2

)T ∣∣∣∣∣ < 1. (5.7)

Let θ = |µ∗2|+ ∆ and assume that |µ∗2| < 2T . If required that
∣∣∣∣ θ

1 + θ
+

1

1 + θ
µ∗2

∣∣∣∣ < 1

2
, which

is equivalent to−1

3
(1 + |µ∗2|) < ∆ < 1 + |µ∗2| , or θ ∈

(
2

3
|µ∗2| −

1

3
, 2 |µ∗2|+ 1

)
. Now, if θ <

2

3
|µ∗2| −

1

3
then 2θ ≤ 2 |µ∗2|+ 1. Therefore, choosing θ = 2k subsequently for k = 1, 2, ...

we are sure that for some k we get θ ∈
(

2

3
|µ∗2| −

1

3
, 2 |µ∗2|+ 1

)
, then the condition (5.7)

will be satisfied. Thus, the grid for sorting parameter θ should be chosen rather coarse. This
justifies the procedure we used in our examples: running (5.6) with small values for θ and then
doubling them until obtaining required cycles. To our surprise, the procedure turns out to be

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



38 D. DMITRISHIN, E. IACOB, A. STOKOLOS

quite efficient. However, this idea has been successfully used recently in adaptive interior-
point methods by Lesaja and Potra [18].

Therefore, the scheme (5.6) allows finding cycles both with small multipliers and with
large ones. In general, the large value of multipliers is not the main obstacle. More challenging
is the problem of small basins of attraction for long cycles. Therefore, it is convenient to either
select a dense grid for initial values or use a sufficiently large number of iterations so that
the point xn would fall into the desired basin of attraction. One can achieve any acceptable
accuracy in determining the cyclic point.

We have experimented with various number of cycles (28, 50, 101, etc.) using Maple and
Python. All the results we include in the subsequent subsections are for T = 101 (sections
5.1–5.13), with a precision of 250 decimals. More results, including the Python code to
replicate the results, can be found in [10].

5.1. Logistic mapping

The logistic mapping
xn+1 = hxn (1− xn) , (5.8)

is, perhaps, the most popular example. Let us consider the case h = 3.99, T = 101. Figure 5.1
illustrates one of its numerous T = 101-cycles.

Fig. 5.1. A 101-cycle of the logistic mapping (5.8).

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



AVERAGE PREDICTIVE CONTROL 39

5.2. Triangular mapping

Our next example is the triangular mapping:

xn+1 = h (1− |2xn − 1|) , h = 0.99. (5.9)

Figure 5.2 shows a T = 101-cyclic point.

Fig. 5.2. A 101-cycle of triangular system (5.9).

5.3. Burgers mapping

For the Burgers mapping:

xn+1 = axn − y2n, yn+1 = byn + xnyn, a = 0.75, b = 1.75 (5.10)

a 101-cyclic point is illustrated in Figure 5.3.

5.4. Tinkerbell mapping

The Tinkerbell mapping:

xn+1 = x2n − y2n + axn + byn, yn+1 = 2xnyn + cxn + dyn, a = 0.9, b = −0.6, c = 2.0, d = 0.5
(5.11)

has a 101-cyclic point illustrated in Figure 5.4.

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



40 D. DMITRISHIN, E. IACOB, A. STOKOLOS

Fig. 5.3. A 101-cycle of the Burgers system (5.10).

Fig. 5.4. A 101-cycle of the Tinkerbell mapping (5.11).

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



AVERAGE PREDICTIVE CONTROL 41

5.5. Gingerbredman mapping

The Gingerbredman mapping:

xn+1 = 1 + |xn| − yn, yn+1 = xn (5.12)

has a 101-cyclic point represented in Figure 5.5.

Fig. 5.5. A 101-cycle of the Gingerbredman mapping (5.12).

5.6. Prey-predator mapping

For the prey-predator mapping:

xn+1 = xn exp (a (1− xn)− byn) , yn+1 = xn (1− exp (−cyn)) , a = 3, b = 5, c = 5
(5.13)

a corresponding 101-cyclic point is illustrated in Figure 5.6.

5.7. Delayed logistic mapping

Figure 5.7 shows a 101-cyclic point of the delayed logistic mapping:

xn+1 = hxn (1− yn) , yn+1 = xn, h = 2.27 (5.14)

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



42 D. DMITRISHIN, E. IACOB, A. STOKOLOS

Fig. 5.6. A 101-cycle of the prey-predator system (5.13).

Fig. 5.7. A 101-cycle of the delayed logistic mapping system (5.14).

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



AVERAGE PREDICTIVE CONTROL 43

5.8. Hénon mapping

The Hénon mapping:

xn+1 = 1 + ax2n + yn, yn+1 = bxn, a = −1.40000001, b = 0.30000002 (5.15)

has a 101-cyclic point represented in Figure 5.8.

Fig. 5.8. A 101-cycle of Hénon system (5.15).

5.9. Elhadj-Sprott mapping

The Elhadj-Sprott mapping:

xn+1 = 1 + a sinxn + byn, yn+1 = xn, a = −4.0, b = 0.9 (5.16)

has a 101-cyclic point illustrated in Figure 5.9.

5.10. Lozi mapping

The Lozi mapping:

xn+1 = 1 + a |xn|+ byn, yn+1 = xn, a = −1.7, b = 0.5 (5.17)

has a 101-cyclic point shown in Figure 5.10.

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



44 D. DMITRISHIN, E. IACOB, A. STOKOLOS

Fig. 5.9. A 101-cycle of the Elhadge-Sprott system (5.16).

Fig. 5.10. A 101-cycle of the Lozi system (5.17).

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



AVERAGE PREDICTIVE CONTROL 45

5.11. Ikeda mapping

The Ikeda mapping is given by the equations:

xn+1 = 1 + u (xn cos τn − yn sin τn) , yn+1 = u (xn sin τn + yn cos τn) , (5.18)

where u = 0.9, τn = 0.4− 6

1 + x2n + y2n
.

The mapping has a 101-cyclic point illustrated in Figure 5.11.

Fig. 5.11. A 101-cycle of the Ikeda system (5.18).

5.12. Holmes cubic mapping

The Holmes cubic mapping:

xn+1 = yn, yn+1 = axn + byn − y3n, a = −0.2, b = 2.77 (5.19)

has a 101-cyclic point shown in Figure 5.12.

5.13. Multihorseshoe mapping

For the Multihorseshoe mapping [14]:

xn+1 = xne
a−0.8xn−0.2yn

yn+1 = yn(0.2xn + 0.8yn)eb−0.2xn−0.8yn

a = 3.0, b = 3.0
(5.20)

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



46 D. DMITRISHIN, E. IACOB, A. STOKOLOS

Fig. 5.12. A 101-cycle of the Holmes cubic system (5.19).

we show a 101-cyclic point in Figure 5.13.

Fig. 5.13. A 101-cycle of the Multihorseshoe system (6.21).

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



AVERAGE PREDICTIVE CONTROL 47

6. CONCLUSION

This article deals with the problem of stabilization for nonlinear systems of two categories:
those unstable and those with a priori unknown periodic orbits at discrete time. A well-known
method of stabilizing controls, called the predictive control method, first proposed by B.T.
Polyak, have been thoroughly investigated in this work. We have found that this method has
several disadvantages: it is necessary to know the cycle exact multiplier or its sufficiently
accurate estimate even in the scalar case; in the vector case, one must know the whole cycle
Jacobi matrix; consequently, the proposed control does not have the required robustness with
respect to the system parameters perturbations; the control gain coefficients have different
signs, which can trigger the initial system multiplier’s shifting beyond the central unit circle
(where it lies) when applying the control; therefore, the gain coefficients must be small, and,
in order to evaluate them at every instance, we need to know the multiplier’s value.

All these shortcomings imply the necessity to modify the predictive control method. We
propose not only to use the first and last iterations of the original mapping, but also all
previous ones, by considering their linear combination. This linear combination’s coefficients
are sought as being of a special polynomial, characterized by certain properties. As a result,
it was possible for us to extend the predictive control scope. In addition, if the coefficients are
non-negative, then for the initial system cycle multipliers lying in the central unit circle the
corresponding multipliers of the control system cycle become closer to zero. An algorithm is
given as a special case of this method application, for finding a cycle of a given length when
its multipliers are known.

One of the possible directions for future research is related to investigating new control
schemes that combine the use of control system previous states and the initial system
predicted states, i.e. the predictive control shall be considered together with the semi-linear
control [8] as follows: 

Xn =
N1∑
j=1

ajxn−jT+T

Yn =
N2∑
j=1

bjxn−jT+1

F (x) =
N3∑
j=1

θjf
((j−1)T+1)(x)

xn+1 = (1− γ)F (Xn) + γYn

(6.21)

where
N1∑
j=1

aj = 1,
N2∑
j=1

bj = 1,
N3∑
j=1

θj = 1. Clearly, the T -cycles of systems (2.1) and (6.21)

coincide. The conditions of the system (6.21) T -cycle local asymptotic stability can be
formulated as

µj [r (µj)]
T ∈

(
C \ Φ

(
D
))∗

, j = 1, . . . ,m,

Φ(z) = (1− γ)T
z(q(z))T

(1− γp(z))T
, q(z) =

N1∑
j=1

ajz
j−1, p(z) =

N2∑
j=1

bjz
j−1,

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



48 D. DMITRISHIN, E. IACOB, A. STOKOLOS

where C is an extended complex plane, and the asterisk denotes the reciprocal operation:

(z)∗ =
1

z
.

The semilinear control method (when N3 = 1 in (6.21)) has also certain disadvantages
[8,9]. Further studies shall aim to eliminating (reducing) inherent disadvantages of predictive
control and semi-linear control, synthesizing these approaches together.

REFERENCES

1. Andrievsky B.R., Fradkov A.L., Control of Chaos: Methods and Applications. I.
Methods, Avtomat. i Telemekh., (2003), no. 5, 3-45.

2. Bartuccelli M., Constantin P., Doering C.R., Gibbon J.D., Gisselfalt M., Hard turbulence
in a finite dimensional dynamical system, Phys. Lett. A 142 (6-7) (1987) 349-356.

3. Chen G., Dong X., From chaos to order: Methodologies, Perspectives and Application.
World Scientific, Singapore (1999)

4. Davidchack R.L., Lai Y.-C., Klebanoff A., Bollt E.M., Towards complete detection of
unstable periodic orbits in chaotic systems, Phys. Lett. A 287 (2001), 99-104.

5. Devaney R.L., An Introduction to Chaotic Dynamical Systems. New York: Addison-
Wesley Publ. Co., Second Edition., (1993)

6. Dmitrishin D., Hagelstein P., Khamitova A., and Stokolos A., Limitations of Robust
Stability of a Linear Delayed Feedback Control, SIAM J. Control Optim. 56, (2018),
148-157.

7. Dmitrishin D. and Khamitova A., Methods of harmonic analysis in nonlinear dynamics,
Comptes Rendus Mathematique, Volume 351, Issues 9-10, 367 - 370 (2013)

8. Dmitrishin D., Khamitova A. and Stokolos A., Fejér polynomials and chaos, Springer
Proceedings in Mathematics and Statistics, 108 (2014), 49–75.

9. Dmitrishin D., Skrinnik I., Lesaja G., Stokolos A., A new method for finding cycles by
semilinear control, Physics Letters A, 383 (2019), 1871 – 1878.

10. Dmitrishin D., Iacob E., Stokolos A., Average predictive control for nonlinear discrete
dynamical systems, arXiv:1906.02925v2.

11. Elaydi S., An Introduction to Difference Equations, ISBN 0-387-23059-9. Springer,
New York (2005), 539p.

12. Galias Z., Rigorous investigations of Ikeda map by means of interval arithmetic.
Nonlinearity, (2002) 15:1759-1779.

13. Jackson E.A., Perspectives of Nonlinear Dinamics. Vol. I, II, - Cambridge Univ. Press,
Cambridge, 1980, 1990 Chaos II, ed. Hao Bai-Lin. - World Sci., (1990).

14. Joshi, Y.A. and Blackmore, D., Strange attractors for asymptotically zero maps. Chaos,
Solitons & Fractals Volume 68, 2014, 123-138.

15. Kittel A, Parisi J, Pyragas K., Delayed feedback-control of chaos by self-adapted delay-
time. Phys Lett A 1995, Vol. 198 (5-6),433-436.

16. Kuntsevich A.V., Kuntsevich V.M., Estimates of Stable Limit Cycles of Nonlinear
Discrete Systems, Journal of Automation and Information Sciences, Volume 44, Issue
9, (2012), 1-10.

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)



AVERAGE PREDICTIVE CONTROL 49

17. Kuntsevich A.V., Kuntsevich V.M., Invariant Sets for Families of Linear and Nonlinear
Discrete Systems with Bounded Disturbances, Automation and Remote Control, Vol.
73, No. 1, (2012), 83-96.

18. Lesaja G., Potra F., Adaptive full newton-step infeasible interior-point method for
sufficient horizontal LCP, Optimization Methods and Software, 21(1), (2018), doi:
10.1080/10556788.2018.1546857

19. Miller J.R., Yorke J.A. Finding all periodic orbits of maps using Newton methods: sizes
of basins. Physica D (2000);135:195-211.

20. Morgul O., Further stability results for a generalization of delayed feedback control,
Nonlinear Dynamics, 1-8 (2012).

21. Morgul O., On the stability of delayed feedback controllers, Phys. Lett. A 314 (2003)
278-285.

22. Ott E., Grebodgi C., Yorke J.A., Controlling chaos. Phys. Rev. Lett. 64, 1196-1199
(1990).

23. Polyak B. T., Stabilizing chaos with predictive control, Automation and Remote
Control, Vol. 66 (11), (2005), 1791-1804.

24. Polyak B.T., Gryazina E.N., Iterations of perturbed tent maps with applications to chaos
control, IFAC Proceedings, Volume 39, Issue 8, 2006, pp 125-129.

25. Pyragas K., Continuous control of chaos by self controlling feedback. Phys. Rev. Lett.
A 170, 421-428 (1992).

26. Sprott, J. C., Chaos and Time-Series Analysis, Oxford University Press, Oxford, Vol. 69
(2003).

27. Shalby L., Predictive feedback control method for stabilization of continuous time
systems. Advances in Systems Science and Applications. 17 (2017), pp 1-13.

28. Ushio T., Limitation of delayed feedback control in nonlinear discrete-time systems.
IEEE Trans Circuits Syst I 1996; Vol. 43 (9), 815-816.

29. Ushio T., Yamamoto, S., Prediction-based control of chaos. Phys. Letts. A. 264(1):439-
446(1999).

30. Vieira de S.M., Lichtenberg A.J. Controlling chaos using nonlinear feedback with delay.
Phys. Rev. E 54, 1200-1207 (1996).

31. Ypma T.J., Historical Development of the Newton-Raphson Method. SIAM Rev.,
(1995) 37:531-551.

32. Zhu J., Tian Y.P., Necessary and sufficient conditions for stabilizability of discrete-time
systems via delayed feedback control. Phys Lett A 2005; Vol. 343 (1-3), 95-107.

Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020)


	Introduction
	Problem statement
	Constructing the Jacobi matrix for a controlled system
	Main result
	Case M={1,…,m}
	Case M={z: Rez0}D
	Case M=D{*}, |*|>1 15
	Case M=D{*,*}, |*|>1
	Case T=1, M="4262304 -*,*"5263305 , M="4262304 -*,1"5263305 
	The general case

	Examples
	Logistic mapping
	Triangular mapping
	Burgers mapping
	Tinkerbell mapping
	Gingerbredman mapping
	Prey-predator mapping
	Delayed logistic mapping
	Hénon mapping
	Elhadj-Sprott mapping
	Lozi mapping
	Ikeda mapping
	Holmes cubic mapping
	Multihorseshoe mapping

	Conclusion