1

Flatness-based Deformation Control of an Euler-Bernoulli Beam with
In-domain Actuation
Amir Badkoubeh, Jun Zheng, and Guchuan Zhu

arXiv:1506.01434v3 [math.OC] 18 Nov 2015

Abstract
This paper addresses the problem of deformation control of an Euler-Bernoulli beam with in-domain actuation. The proposed
control scheme consists in first relating the system model described by an inhomogeneous partial differential equation to a target
system under a standard boundary control form. Then, a combination of closed-loop feedback control and flatness-based motion
planning is used for stabilizing the closed-loop system around reference trajectories. The validity of the proposed method is
assessed through well-posedness and stability analysis of the considered systems. The performance of the developed control
scheme is demonstrated through numerical simulations of a representative micro-beam.

I. I NTRODUCTION
The present work addresses the control of an Euler-Bernoulli beam with in-domain actuation described by an inhomogeneous
partial differential equation (PDE). A motivating example of such a problem arises from the deformation control of microbeams, shown in Fig. 1. This device is a simplified case of deformable micro-mirrors that are extensively used in adaptive optics
[1], [2]. Due to technological restrictions in the design, the fabrication and the operation of micro-devices, design methods that
may lead to control structures with a large number of sensors and actuators are not applicable to microsystems with currently
available technologies.
One of the standard methods to deal with the control of inhomogeneous PDEs is to discretize the PDE model in space to
obtain a system of lumped ordinary differential equations (ODEs) [3], [4]. Then, a variety of techniques developed for the
control of finite-dimensional systems can be applied. However, in addition to the possible instability due to the phenomenon
of spillover [5], the increase of modeling accuracy may lead to high-dimensional and complex feedback control structures,
requiring a considerable number of actuators and sensors for the implementation. Therefore, it is of great interest to directly
deal with the control of PDE models, which may result in control schemes with simple structures.
There exists an extensive literature on the control of flexible beams described by PDEs in which the majority of the reported
work deal with the stabilization problem by means of boundary control (see, e.g., [3], [6], [7], [8], [9], [10]). Placing actuators
in the domain of the system will lead to inhomogeneous PDEs [3], [11], [12]. The stabilization problem of Euler-Bernoulli
beams with in-domain actuation is considered in, e.g., [11], [13], [14], [15], [16]. The deformation control of beams has gained
in recent years an increasing attention and impotence (see, e.g., [17], [18], [19]).
In the present work, we employ the method of flatness-based control, which has been applied to a variety of infinitedimensional systems (see, e.g., [18], [20], [21], [22], [23], [24], [25], [26], [27]). Nevertheless, applying this tool to systems
controlled by multiple in-domain actuators leads essentially to a multiple-input multiple-output (MIMO) problem, which still
remains a challenging topic. It is noticed that a scheme proposed recently in [18], [19] tackles this issue by utilizing the
Weierstrass-factorized representation of the spectrum of the input-output dynamics. Nevertheless, an early truncation is still
required in order to obtain a finite-dimensional input-output map, which is essential to apply the technique of flat systems.
The control scheme developed in this paper is based on an approach aimed at avoiding early truncations in control
design procedure. Specifically, it consists in first relating the original inhomogeneous model with pointwise actuation to a
system in a standard boundary control form to which the technique of flatness-based motion planning may be applied for
feedforward control. By using the technique of lifting, this boundary-controlled system can be transformed into a regularized
inhomogeneous system driven by sufficiently smooth functions. It is shown that the steady-state solution of the regularized
system can approximate that of the original system. A standard closed-loop feedback control is used to stabilize the original
inhomogeneous system. Moreover, the MIMO issue is solved by splitting the reference trajectory into a set of sub-trajectories
based on an essential property of the Green’s function. Well-posedness and stability analysis of the considered systems is also
carried out, which is essential for the validation of the developed control scheme. It is worth noting that the approach proposed
in this paper does not involve any early truncation. Moreover, it is shown that the invertibility of the resulting input-output
map can be guaranteed a priori. This allows for the control design to be carried out directly with the original PDE model in
a systematic manner.
This paper is a preprint of a paper submitted to IET Control Theory and Applications and is subject to Institution of Engineering and Technology Copyright.
If accepted, the copy of record will be available at IET Digital Library
A. Badkoubeh with the the Department of Electrical Engineering, Ecole Polytechnique de Montreal, P.O. Box 6079, Station Centre-Ville, Montreal, QC,
Canada H3T 1J4. E-mail: amir.badkoubeh@polymtl.ca
J. Zheng is with the Department of Basic Course, Southwest Jiaotong University, Emeishan, Sichuan 614202, China. E-mail: zheng@xx
G. Zhu is with the Department of Electrical Engineering, Ecole Polytechnique Montreal, P.O. Box 6079, Station Centre-Ville, Montreal, QC, Canada H3T
1J4. E-mail: guchuan.zhu@polymtl.ca

2

Fig. 1.

Schematic of the deformable microbeam.

The remainder of the paper is organized as follows. Section II describes the model of a representative in-domain actuated
deformable beam. Section III introduces the relationship between the interior actuation and the boundary control and examines
the well-posedness of the considered systems. Section IV addresses the feedback control and stability issues. Section V deals
with motion planning for feedforward control. A simulation study is carried out in Section VI, and, finally, some concluding
remarks are presented in Section VII.
II. M ODEL OF THE I N - DOMAIN ACTUATED E ULER -B ERNOULLI B EAM
As shown in Fig. 1, the considered structure consists of a continuous flexible beam and an array of actuators connected to
the beam via rigid spots. As the dimension of the spots connecting the actuators to the continuous surface is much smaller
than the extent of the beam, the effect of the force generated by actuators can be considered as a pointwise control represented
by Dirac delta functions concentrating at rigid spots. Furthermore, it is supposed that the structure is supported at the left
end, terminated to the frame, and is suspended at the right end, vertically suspended on a moving actuator, corresponding to
a simply supported-shear hinged configuration [15].
As the developed scheme uses one actuator for stabilizing feedback control, for notational convenience we consider a setup
with N + 1 actuators. The displacement of the beam at the position x and the instance t is denoted by w(x, t). The derivatives
of w with respect to its variables are denoted by wx and wt , respectively. Suppose that the mass density and the flexural
rigidity of the beam are constant. Then, the dynamic transversal displacement of the beam with pointwise actuators located at
{x1 , x2 , · · · , xN , xN +1 } in a normalized coordinate can be described by the following PDE [14], [15], [28]:
wtt (x, t) + wxxxx (x, t) =

N
+1
X

αj (t)δ(x − xj ), x ∈ (0, 1) , t > 0,

(1a)

j=1

w(0, t) = wx (1, t) = wxx (0, t) = wxxx (1, t) = 0, t > 0,

(1b)

w(x, 0) = h0 (x), wt (x, 0) = h1 (x), x ∈ (0, 1),

(1c)

where x is a normalized variable spanned over the domain (0, 1), δ(x − xj ) is the Dirac mass concentrated at the point
xj ∈ (0, 1] and αj : t 7→ R, j = 1, . . . , N + 1, are the control signals. h0 and h1 represent the initial value of the beam.
Without loss of generality, we assume that 0 < x1 < x2 < · · · < xN < xN +1 ≤ 1.
Remark 1: The considered problem can also be modeled as a serially connected beam [15]. Indeed, these two models are
equivalent in the sense that they lead to the same abstract linear system in variational form (see, e.g., [29]).
The control objective is to steer the beam to a desired form in steady state via in-domain pointwise actuation.
The model given in (1) is a nonstandard PDE due to the unbounded inputs, represented by Dirac functions, on the right-hand
side of (1a). Therefore, we need to invoke the weak derivative in the theory of distributions [30] and discuss solutions of the
considered problem in a weak sense. To this end, we start by defining the convergence on C0∞ (Ω).
Definition 1: Let Ω be a domain in R. A sequence {ϕj } of functions belonging to C0∞ (Ω) converge to ϕ ∈ C0∞ (Ω), if
(i) there exists K ⊂ Ω such that supp(ϕj − ϕ) ⊂ K for every j, and
(ii) limj→∞ ∂ n ϕj (x) = ∂ n ϕ(x) uniformly on K for all n ≥ 0.
The linear space C0∞ (Ω) having the above property of convergence is called fundamental space, denoted by D(Ω). The
space of all linear continuous functionals on D(Ω), denoted by D 0 (Ω), is called the space of (Schwartz) distributions on Ω,
which is the dual of D(Ω) (see, e.g., Chapter 1 of [31] for more properties of D(Ω) and D 0 (Ω)).
Denote the set Φ by

Φ = φ ∈ H 2 (0, 1) ; φ(0) = φx (1) = 0 .
(2)
We give the definition of weak solution of (1).
Definition 2: Let T > 0 and αi ∈ L2 (0, T ) for i = 1, 2, ..., N + 1. Let h0 ∈ Φ, h1 ∈ L2 (0, 1). A weak solution to the
problem (1) is a function w ∈ C([0, T ] ; Φ) ∩ C 1 ([0, T ] ; L2 (0, 1)), satisfying
w(x, 0) = h0 (x),

x ∈ (0, 1) ,

3

such that, for every v ∈ C 1 ([0, T ] ; D(0, 1)), one has for almost every τ ∈ [0, T ]
Z 1
Z 1
wt (x, τ )v(x, τ )dx −
h1 (x)v(x, 0)dx
0
0
Z τZ 1
wt (x, t)vt (x, t)dxdt
−
0
0
Z τZ 1
wxx (x, t)vxx (x, t)dxdt
+
=

0
0
N
+1 Z τ Z 1
X
i=1

0

(3)

αi (t)δ(x − xi )v(x, t)dxdt.

0

III. I N - DOMAIN ACTUATION D ESIGN VIA B OUNDARY C ONTROL
A. Relating In-domain Actuation to Boundary Control
In the scheme developed in the present work, we use N feedforward control signals to deform the beam, while dedicating
the feedback stabilizing control to one actuator. For notational simplicity, we assign the feedback control to the (N + 1)th
actuator. Hence, the dynamic of the desired trajectory driven by the feedforward control can be expressed as
d
d
wtt
(x, t) + wxxxx
(x, t) =

N
X

αj (t)δ(x − xj ), x ∈ (0, 1) , t > 0,

(4a)

j=1
d
d
wd (0, t) = wxd (1, t) = wxx
(0, t) = wxxx
(1, t) = 0, t > 0,
d

w (x, 0) =

wtd (x, 0)

= 0, x ∈ (0, 1) .

(4b)
(4c)

A weak solution of (4) can be defined in a similar way as that of (1) described in Definition 2. Note that it is shown later that
the initial condition given in (1c) will be captured by the regulation error dynamics. Therefore, feedforward control design can
be carried out based on System (4) with zero-initial conditions.
Due to the fact that the model given in R(4) is driven by unbounded inputs, we will apply a sequence of cutoff functions
+∞
ϕ(x − xj ) having unit mass and satisfying −∞ (x − xj )ϕ(x − xj )dx = 0, also called blobs [32], to approximate δ(x − xj ) in
the sense of distributions. We will then show that in steady state, wd can be approximated by a sufficiently smooth function.
Specifically, we consider first the following boundary controlled PDE:
utt (x, t) + uxxxx (x, t) = 0, x ∈ (0, 1) , t > 0,

(5a)

u(0, t) = ux (1, t) = uxx (0, t) = 0, t > 0,

(5b)

uxxx (1, t) = g(t), t > 0,

(5c)

u(x, 0) = ut (x, 0) = 0, x ∈ (0, 1) ,

(5d)

PN

3
where g(t) =
j=1 gj (t). Throughout this paper, without special statements, we assume that gj (t) ∈ C ([0, +∞)) and
gj (0) = ġj (0) = 0 for j = 1, 2, ..., N . Notice that the motivation behind considering the system of the form given in (5) is
that it allows employing the techniques of boundary control for feedforward control design while avoiding early truncations
of dynamic model and/or controller structure.
A weak solution to (5) can be defined in a similar way as that of (1) described in Definition 2.
Definition 3: Let T > 0. A weak solution to the problem (5) is a function u ∈ C([0, T ] ; Φ)∩C 1 ([0, T ] ; L2 (0, 1)) satisfying

u(x, 0) = ut (x, 0) = 0,

x ∈ (0, 1) ,

1

such that, for every v ∈ C ([0, T ] ; Φ), one has for almost every τ ∈ [0, T ]
Z 1
Z τZ 1
ut (x, τ )v(x, τ )dx −
ut (x, t)vt (x, t)dxdt
0
0
0
Z τ
Z τZ 1
+
g(t)v(1, t)dt +
uxx (x, t)vxx (x, t)dxdt = 0.
0

Let ψ(x, t) = u(x, t) −

PN

j=1 gj (t)Hj (x),

0

(6)

0

where Hj (x), j = 1, 2, ..., N , are defined in [0, 1] satisfying

Hjxxxx (x) = ϕ(x − xj ),

(7a)

Hj (0) = Hjx (1) = Hjxx (0) = 0, Hjxxx (1) = 1,

(7b)

4

with ϕ(x) ∈ L2 (R). Then by lifting, (5) can be transformed into the following PDE with zero boundary conditions
ψtt (x, t) + ψxxxx (x, t)
=−

N
X

g̈j (t)Hj (x) −

j=1

N
X

gj (t)Hjxxxx (x), x ∈ (0, 1) , t > 0,

(8a)

j=1

ψ(0, t) = ψx (1, t) = ψxx (0, t) = ψxxx (1, t) = 0,

(8b)

ψ(x, 0) = ψt (x, 0) = 0, x ∈ (0, 1) .

(8c)

In order to relate (4) to (5), we first establish a relationship between (4) and (8), especially in steady state. Let αj (t) = −gj (t)
and suppose that limt→∞ αj (t) = αj and limt→∞ gj (t) = g j for j = 1, . . . , N . We have then in steady state:
ψ xxxx (x) = −

N
X

g j Hjxxxx (x)

j=1

=

N
X

αj ϕ(x − xj ), x ∈ (0, 1) , t > 0,

(9a)

j=1

ψ(0) = ψ x (1) = ψ xx (0) = ψ xxx (1) = 0.

(9b)

For j = 1, 2, ..., N , given a sequence of blobs {ϕm (x − xj )}, we seek a sequence of functions {Hjm (x)} such that
m
Hjxxxx
(x) = ϕm (x − xj )

(10)

with Hjm (x) satisfying (7b) and ϕm (x − xj ) → δ(x − xj ) in D 0 (0, 1) as m → +∞. Therefore, considering (9) and the
steady-state model of (4), we have ψ m → wd in C 1 ([0, 1]) as m → +∞ (see Theorem 2). Hence, for feedforward control
design, we may consider the systems (8) and (9).
Lemma 1: [32] Let ϕm (x) ∈ L2 (R) be defined by
ϕm (x) =

1 sin mx
.
π x

(11)

Then ϕm (x) has the following properties:
Z +∞
Z +∞
1 sin mx
dx = 1.
(i)
ϕm (x)dx =
x
−∞
−∞ π
0
(ii) ϕm (x) → δ(x) in D (R) as m → ∞.
By taking ϕm given in (11) as the input to (10), we get
1
1
Hjm (x) = x3 − x
6 Z Z2 Z
x
z
+
0

1

y

Z

t

ϕm (s − xj )dsdtdydz.

0

1

Theorem 2: Let ϕm (x) be defined as in Lemma 1. Assume that αj (t) = −gj (t) tends to αj = −g j as t → ∞ for all
m
m
j = 1, 2, ..., N . Denote by ψ j and by wdj the steady-state solutions of System (4) and System (8), respectively. Then ψ j → wdj
in C 1 ([0, 1]) as m → +∞ for j = 1, 2, ..., N .
Proof: In the steady state, we have
m

m
ψ jxxxx (x) = −g j Hjxxxx
(x) = αj ϕm (x − xj ),

(12a)

m
ψ j (0)

(12b)

=

m
ψ jx (1)

=

m
ψ jxx (0)

=

m
ψ jxxx (1)

=

wdjxxx (1)

= 0,

and
wdjxxxx (x) = αj δ(x − xj ),
wdj (0)

=

wdjx (1)

=

wdjxx (0)

(13a)
= 0.

Taking v(x) ∈ D(0, 1), with v(0) = vx (1) = vxx (0) = 0, as a test function and integrating by parts, we get
Z 1

m
ψ jx (x) − wdjx (x) vxxx (x)dx
0
Z 1
=αj
(ϕm (x − xj ) − δ(x − xj )) v(x)dx.
0

(13b)

5
m

Since ϕm (x − xj ) → δ(x − xj ) in the sense of distributions as m → +∞ and vxxx ∈ D(0, 1), it follows that ψ jx → wdjx
m
in the sense of distributions as m → +∞ for j = 1, 2, ..., N . Furthermore, as ψ jx ∈ L1 (0, 1) and wdjx ∈ L1 (0, 1), we have
m
m
ψ jx → wdjx a.e. in (0, 1) (see Lemma 3.31 of [31], page 74). Then by the continuity of ψ jx and wdjx , we conclude that
m
m
ψ jx → wdjx pointwise in (0, 1). Therefore ψ j → wdj in C 1 ([0, 1]) as m → +∞ for j = 1, 2, ..., N .
Remark 2: Note that describing the control acting on discrete points in the domain by Dirac delta functions is a mathematical
abstraction, which allows for the beam equation to possess some fundamental properties, in particular the exact controllability
(see, e.g., [33]), and facilitates control design. However, the real physical devices cannot generate unbounded control concentrated on a point. From this viewpoint, Theorem 2 provides an assessment of the gap between the performance achievable by
utilizing real physical devices and that predicted by the corresponding idealized mathematical model.
B. Well-posedness of Cauchy Problems
Well-posedness analysis is essential to the approach developed in this work. In this subsection, we establish the existence
and the uniqueness of weak solutions of equations (1), (4), (5) and (8).
Theorem 3: The following statements hold true:
(i) Assume αj ∈ L2 (0, T ) for j = 1, 2, ..., N + 1. Let h0 ∈ Φ, h1 ∈ L2 (0, 1) and T > 0. Then System (1) and System
(4) has a unique weak solution w ∈ C([0, T ] ; Φ) ∩ C 1 ([0, T ] ; L2 (0, 1)) and wd ∈ C([0, T ] ; Φ) ∩ C 1 ([0, T ] ; L2 (0, 1)),
respectively.
(ii) Let T > 0 and g ∈ C 2 ([0, T ]). Then System (5) has a unique weak solution u ∈ C([0, T ] ; Φ) ∩ C 1 ([0, T ] ; L2 (0, 1)).
Furthermore, if g ∈ C 3 ([0, T ]) and Hjm (x), j = 1, 2, ..., N , is defined as in (7b), then System (8) has a unique solution
ψm ∈ C([0, T ] ; Φ) ∩ C 1 ([0, T ] ; L2 (0, 1)).
Proof: The proof of (i) can be proceeded step by step as in Proposition 3.1 of [11]. We prove the first result of (ii) and
the second part can be proceeded in the same way. We assume first g(t) ∈ C 3 ([0, T ]) and consider the following system:


1
1
x − x3 g̈(t), x ∈ Ω, t > 0,
(14a)
νtt (x, t) + νxxxx (x, t) =
2
6
ν(0, t) = νx (1, t) = νxx (0, t) = νxxx (1, t) = 0,
(14b)
ν(x, 0) = 0, νt (x, 0) = 0, x ∈ Ω.

(14c)

4
Let X = Φ × L2 (0, 1) and H(0)
(0, 1) = {u ∈ H 4 (0, 1); u(0) = ux (1) = uxx (0) = uxxx (1) = 0}. Define the inner product
R1
on X by h(u1 , v1 ), (u2 , v2 )iX = 0 (u1xx u2xx + v1 v2 )dx. Define the subspace D(A0 ) ⊂ X by D(A0 ) = {(u, v); (u, v) ∈
4
H(0)
(0, 1) × Φ}, with the corresponding operator A0 : D(A0 ) → X defined as:
  

u
v
A0
=
.
v
−uxxxx

One may easily check that D(A0 ) is dense in X, A0 is closed, A∗0 = −A0 , and hA0 z, ziX = 0. Thus, by Stone’s theorem, A0
generates a semigroup of isometries on X. Based on a classical result on perturbations of linear evolution equations (see, e.g.,
Theorem 1.5, Chapter 6, page 187, [34]) there exists uniquely z = (z1 , z2 ) ∈ C 1 ([0, T ] ; X) ∩ C([0, T ] ; D(A0 )), such that


0 

dz
 , z(·, 0) = (0, 0),
1
= A0 z +  1
x − x3 g̈(t)
dt
2
6
4
which implies that (14) has a unique solution ν = z1 ∈ C([0, T ] ; H(0)
(0, 1)) ∩ C 1 ([0, T ] ; Φ) in the usual sense. Particularly,

1
2
ν ∈ C([0, T ] ; Φ) ∩ C ([0, T ] ; L (0, 1)) is a weak solution. A direct computation shows that u = ν − 21 x − 16 x3 g(t) is a
solution of (5) in the usual sense and, in particular, it is a weak solution.

Now for g ∈ C 1 ([0, T ]), let gn ∈ C 3 ([0, T ]) such that gn → g in C 1 ([0, T ]). Consider un = νn − 12 x − 61 x3 gn (t), where
νn is the solution of (14) corresponding to the data gn (t). Then arguing as above and taking limit as in Chapter 2 of [35], we
can obtain the existence and the uniqueness of a weak solution of (5).

IV. F EEDBACK C ONTROL AND S TABILITY OF THE I NHOMOGENEOUS S YSTEM
The validity of the proposed scheme requires a suitable closed-loop control that can guarantee the stability of the original
inhomogeneous system. As the (N + 1)th actuator is dedicated to stabilizing control, (1a) can be written as
wtt (x, t) + wxxxx (x, t) − αN +1 (t)δ(x − xN +1 )
=

N
X
j=1

αj (t)δ(x − xj ), x ∈ (0, 1), t > 0.

(15)

6

Suppose further that the feedback control is taken as [11], [15]
αN +1 (t) = −kwt (xN +1 , t) ,

(16)

where k is a positive-valued constant. Then in closed-loop, (15) becomes
wtt (x, t) + wxxxx (x, t) + kwt (x, t) δ(x − xN +1 )
=

N
X

αj (t)δ(x − xj ), x ∈ (0, 1), t > 0.

(17)

j=1

Let D(A) = {(w, v); (w, v) ∈ (H 2 (0, 1) ∩ (H 4 (0, x1 ) ∪ H 4 (x1 , x2 ) ∪ · · · ∪ H 4 (xN , xN +1 ) ∪ H 4 (xN +1 , 1))) × H 2 (0, 1), w(0) =
wx (1) = wxx (0) = wxxx (1) = 0, v(0) = vx (1) = 0} and X be defined as in the proof of Theorem 3. Let T > 0 and
αj ∈ L2 (0, T ) for j = 1, ..., N . Assume (h0 , h1 ) ∈ X. To address the stability of System (17) with the boundary conditions
(1b) and initial conditions (1c), we consider the corresponding linear control system under the following abstract form:
ż = Az + Bα, t > 0,
0

T

z(0) = z = (h0 , h1 ) ,
T

(18a)
(18b)

T

where z = (w, v) , α = (α1 , · · · , αN ) , A: D(A) → X is defined as:
  

w
v
A
=
,
v
−wxxxx − kvδ(x − xN +1 )
with v = wt and B: RN → D0 (A∗ ), where A∗ is the adjoint of A, which is defined as:

T
0
(Bα)y = PN
y, ∀ y ∈ D(A∗ ).
j=1 αj δ(x − xj )

(19)

(20)

Note that B ∗ : D(A∗ ) → RN is defined by
B ∗ y = (y2 (x1 ), . . . , y2 (xN ))T , ∀y =

 
y1
∈ D(A∗ ).
y2

(21)

Let U = RN . Then the solution of the Cauchy problem (18) can be defined as follows (see [35], Definition 2.36, p53):
hz(τ ), y τ iX
=hz 0 , S ∗ (τ )y τ iX +

Z

τ

hα(t), B ∗ S ∗ (τ − t)y τ iU dt,

(22)

0

∀ τ ∈ [0, T ] , ∀ y τ ∈ X.
One may verify, as in Chapter 2 of [35], that the solution defined by (22) is also a weak solution of (17) under the form given
in Definition 2.
Theorem 4: Assume αj ∈ L∞ (0, +∞), j = 1, . . . , N . For any xN +1 ∈ (0, 1] being a rational number with coprime
factorization, there exist positive constants C1 , C2 and λ such that for any 0 ≤ τ ≤ T < ∞, there holds:
kz(τ )kX ≤ C1 e−λτ kz 0 kX + C2 kαkL∞ (0,T ) .

(23)

Proof: The proof is a slight modification of Theorem 2.37 in [35]. First, the admissible property of B can be obtained as
(3.38) in [11]. A direct computation gives
    
 
u
u
u
2
A
,
= −k|v(xN +1 )| ,
∀
∈ D(A).
v
v
v
X
Therefore, A is a dissipative operator. Moreover, it can be directly verified that A is onto and hence, according to Theorems
4.3 and 4.6 of [30] (p. 14-15), it generates a C0 -semigroup of linear contractions S(t) acting on X. Furthermore, S(t) is
exponentially stable if xN +1 ∈ (0, 1) is a rational number with coprime factorization, in particular for xN +1 = 1 (see, e.g.,
[11], [15]), i.e. there exist two positive constants C1 and λ, such that
kS(t)kL(X;X) ≤ C1 e−λt , ∀ t ≥ 0.

7

Then for a weak solution defined in (22), we have that for all 0 ≤ τ ≤ T < ∞:
hz(τ ), y τ iX
0

Z

τ

=hS(τ )z , y iX +

τ

hS(τ − t)Bα(t), y τ iU dt

0
Z

=hS(τ )z 0 , y τ iX +

τ

S(τ − t)Bα(t)dt, y τ

0

≤C1 e

−λτ

0


U

τ

τ

Z

kz kX ky kX + C1 kαkL∞ (0,T ) ky kL∞ (0,T )

τ

e−λ(τ −t) dt

0

C1
kαkL∞ (0,T ) ky τ kL∞
λ
C1
kαkL∞ (0,T ) ky τ kH 1 .
≤C1 e−λτ kz 0 kX ky τ kX +
λ
Applying the interpolation formula
kux kL2 ≤ c (kukL2 + kuxx kL2 ) ,
≤C1 e−λτ kz 0 kX ky τ kX +

where c is a positive constant, and the Sobolev embedding theorem yields
ky τ kH 1 ≤ cky τ kH 2 ≤ cky τ kX .
We have then
τ

hz(τ ), y iX



C1 c
−λτ
0
≤ C1 e
kz kX +
kαkL∞ (0,T ) ky τ kX ,
λ

which implies that (23) holds with C2 = C1 c/λ.
Now consider the regulation error defined as e(x, t) = w(x, t)−wd (x). Denoting ∆αj (x, t) = αj (x, t)−αj (x), j = 1, . . . , N ,
then from (17), (1b), (1c) and the steady-state model of (4), the regulation error dynamics satisfy
ett (x, t) + exxxx (x, t) + ket (x, t) δ(x − xN +1 )
=

N
X

∆αj (t)δ(x − xj ), x ∈ (0, 1), t > 0,

(24a)

j=1

e(0, t) = ex (1, t) = exx (0, t) = exxx (1, t) = 0, t > 0,

(24b)

d

e(x, 0) = e0 (x) = h0 (x) − w (x),
et (x, 0) = e1 (x) = h1 (x),

x ∈ (0, 1).

(24c)

Obviously, the regulation error dynamics are in an identical form as (17) with the same type of boundary conditions. We
can then consider the solution of System (24) defined in the same form given by (22) with ze = (e, et ) and ze0 = (e0 , e1 ).
Corollary 5: Assume that all the conditions in Theorem 4 are fulfilled and ze0 ∈ X. Then there exist positive constants C1 ,
C2 and λ, independent of t, such that for any t ≥ 0, there holds:
kze (t)kX ≤ C1 e−λt kze0 kX + C2 k∆αkL∞ (0,+∞) .

(25)

Moreover, if limt→∞ ∆αj = 0 for all j = 1, . . . , N , then limt→∞ e(x, t) = 0, ∀x ∈ (0, 1).
The feedforward control satisfying the conditions for closed-loop stability and regulation error convergence can be obtained
through motion planning, as presented in the next section.
V. M OTION P LANNING AND F EEDFORWARD C ONTROL
According to the principle of superposition for linear systems, we consider in feedforward control design the dynamics of
System (5) corresponding to the input gj (t) described in the following form:
ujtt (x, t) + ujxxxx (x, t) = 0,

(26a)

uj (0, t) = ujx (1, t) = ujxx (0, t) = 0,

(26b)

ujxxx (1, t) = gj (t),

(26c)

uj (x, 0) = ujt (x, 0) = 0.
The required control signal gj (t) should be designed so that the output of System (26) follows a prescribed function
while guaranteeing that the conditions in Theorem 4 and Corollary 5 are fulfilled.

(26d)
udj (xj , t),

8

Motion planning amounts then to taking udj (xj , t) as the desired output and to generating the full-state trajectory of the
subsystem uj (x, t). By applying a standard Laplace transform-based procedure (see, e.g., [18], [20], [22], [27]), we can obtain
the full-state trajectory with zero initial values expressed in terms of the so-called flat output, yj (t), and its time-derivatives:


1
1
uj (x, t) =
x − x3 yj (t)
2
6
∞
n
X
X
x4k+1
+
(27)
(4k + 1)!(4(n−k) + 2)!
n=1 k=0
!
n
X
x4k+3
(2n)
−
(−1)n yj (t).
(4k + 3)!(4(n−k))!
k=0

Now let yj (t) = y j φj (t), where φj (t) is a smooth function evolving from 0 to 1. The full-state trajectory can then be written
as

uj (x, t) =y j

1
x−
2
∞
X


1 3
x φj (t)
6
n
X

x4k+1
(4k + 1)!(4(n−k) + 2)!
n=1 k=0
!
n
X
x4k+3
(2n)
−
(−1)n φj (t).
(4k + 3)!(4(n−k))!

+ yj

(28)

k=0

The corresponding input can be computed from (26c), which yields
gj (t) = − y j φj (t) + y j

∞
n
X
X
n=1

−

n
X
k=0

k=1

1
(4k)!(4(n−k))!

1
(4k − 2)!(4(n−k) + 2)!
(29)

!
(2n)
(−1)n φj (t).

For set-point control, we need an appropriate class of trajectories enabling a rest-to-rest evolution of the system. A convenient
choice of φj (t) is the following smooth function:

0,
if t ≤ 0

Z t




ε

exp(−1/(τ (1 − τ ))) dτ

0
,
if t ∈ (0, T )
φj (t) =
(30)
Z T


ε

exp(−1/(τ
(1
−
τ
)))
dτ



 0
1,
if t ≥ T
which is known as Gevrey function of order σ = 1 + 1/ε, ε > 0 (see, e.g., [27]).
For the convergence of (28) and (29), we have
Proposition 6: If φj (t) in the basic output y j (t) = y j φj (t) is chosen as a Gevrey function of order 1 < σ < 2, then the
infinite series (28) and (29) are convergent.
The proof of Proposition 6 is based on the bounds of Gevrey functions of order σ given by (see, e.g., [21], [36]):
(k!)σ
(k+1)
∃K, M > 0, ∀k ∈ Z≥0 , ∀t ∈ [t0 , T ] , φj
(t) ≤ M
.
(31)
Kk
The detailed convergence analysis is given in Appendix IX-A.
Remark 3: In general, the Gevrey bounds are unknown, but can be estimated following the way presented in [21]. Furthermore, a symmetric function in the transient phase can be considered to improve convergency analysis [21].
Now let
1
1
Pj (x) = x − x3 ,
2Z Z 6 Z Z
x

z

y

t

Ij,m (x) =
ϕm (s − xj )dsdtdydz,
0
1
0
1
X

n
x4k+1
x4k+3
Φj,n (x) =
−
(−1)n ,
(4k + 1)!(4(n−k) + 2)! (4k + 3)!(4(n−k))!
k=0
X

n
n
X
1
1
Ψj,n (1) =
−
(−1)n .
(4k − 2)!(4(n−k) + 2)!
(4k)!(4(n−k))!
k=1

k=0

9

Set

PN

j=1

ψjm (x, t) = ψm (x, t), m > 0. By the definition of ψ(x, t), we obtain
ψjm (x, t) =y j Ij,m (x)φj (t) + y j

∞
X

(Φj,n (x)

(33)

n=1

+Ψj,n (1)Pj (x) −

(2n)
Ψj,n (1)Ij,m (x)) φj (t).

To compute Ij,m (x), we note that
ϕm (x − xj ) =

∞
1 X (−1)k+1 m2k−1
(x − xj )2k−2 .
π
(2k − 1)!
k=1

Therefore,
Ij,m (x) =


∞
1 X (−1)k+1 m2k−1 (x − xj )2k+2
π
(2k − 1)
(2k + 2)!
k=1

2
x2k
(1 − xj )2k−1 x3
(1 − xj )2k+1 x
j x
−
−
6((2k − 1)!)
2(2k)!
(2k + 1)!
2k+2 
2k−1
x
x2k
x
(1 − xj )
x
j
j
+
+
−
.
(2k)!
2(2k − 1)!
(2k + 2)!

−

(34)

Claim 7: Ij,m (x) given in (34) is convergent for all x, xj ∈ (0, 1) with respect to any fixed m.
The proof of this claim is given in Appendix IX-B.
Based on Claim 7, the series in the right hand side of (33) are convergent and ψjm (x, t) can be expanded by (33) and (34).
m
Moreover, for gj (t) given in (29), ψjm (x, t) tends to ψ j (x) = y j Ij,m (x) as t → ∞. Note that for fixed x, the radius of
convergence of ψ j,m on m is ∞. Thus, we can let m → +∞.
To complete the control design, we need to determine the amplitude of flat outputs y j , j = 1, . . . , N , from the desired shape
w
ed (x) that may not necessarily be a solution of the steady-state beam equation (13). To this end, we propose to interpolate
w
ed (x) using the Green’s function of (13), denoted by G(x, ξ), which is the impulse response corresponding to the input
δ(x − ξ) and is given by



3

 − x + xξj 1 − ξj , 0 ≤ x < ξj ;

6
2
G(x, ξj ) =
3


ξ

x

 − j + ξj x 1 −
,
ξj ≤ x ≤ 1.
6
2
Due to the principle of superposition for linear systems, the solution to (13), wd (x), can be expressed as
Z 1X
N
wd (x) =
G(x, ξ)αj δ(x − ξj )dξ
0 j=1

=

N
X

G(x, ξj )αj .

j=1

Now taking N points on w
ed (x) and letting wd (xj ) = w
ed (xi = j), j = 1, . . . , N , yield


   d
G(x1 , ξ1 ) . . . G(xN , ξ1 )
α1
w
e (x1 )


  ..  
..
..
..
..

,
 .  = 
.
.
.
.
G(x1 , ξN ) . . .

G(xN , ξN )

αN

w
e (xN )

which represents a steady-state input to output map.
Claim 8: The map of (35) is invertible for all xj , ξj ∈ (0, 1), j = 1, . . . , N and xi 6= xj , ξi 6= ξj , if i 6= j.
The proof of this claim is given in Appendix IX-C.
As αj = −g j and limt→∞ gj (t) = g j = −y j for all j = 1, . . . , N , we obtain from Claim 8 that
−1 d
  

y1
G(x1 , ξ1 ) . . . G(xN , ξ1 )
w
e (x1 )
 
 ..  

..
..
..
..
 . =
 
.
.
.
.
.
yN

G(x1 , ξN ) . . .

G(xN , ξN )

(35)

d

d

w
e (xN )

(36)

10

VI. S IMULATION S TUDY
In the simulation study, we consider the deformation control in which the desired shape is given by

2
2
w
ed (x) = − 10−3 e−100(x−0.4) + 2e−100(x−0.6)

2
+3e−400(x−0.7) , x ∈ (0, 1),

(37)

as shown in Fig. 2(a). Note that the maximum amplitude of the shape given by (37) is 3.8 × 10−3 , which represents a typical
micro-structure for which the beam length is of few centimeters and the displacement is of micrometer order.
In order to obtain an exponential closed-loop convergency, the actuator for feedback stabilization is located at the position
xN +1 = 1. To evaluate the effect of the number of actuators to interpolation accuracy, measured by kw
ed (x) − wd (x)kL1 (0,1)
and control effort, we considered 3 setups with, respectively, 8, 12 and 16 actuators evenly distributed in the domain. It can
be seen from Fig. 2 that the setup with 8 actuators exhibits an important interpolation error and the one with 16 actuators
requires a high control effort in spite of a high interpolation accuracy. The setup with 12 actuators provides an appropriate
trade-off between the interpolation accuracy and the required control effort, which is used in control algorithm validation.
A MATLAB Toolbox for dynamic Euler-Bernoulli beams simulation provided in Chapter 14 of [37] is used in numerical
implementation. With this Toolbox, the simulation accuracy can be adjusted by choosing the number of modes used in
2
implementation. In the simulation, we implement System (1) with initial conditions h0 (x, 0) = −3 × 10−3 e−400(x−0.8)
and h1 (x, 0) = 0. As an undamped beam is unstable in open-loop, the controller tuning is started by determining a suitable
value of the closed-loop control gain k. The basic outputs φj (t) used in the simulation are Gevrey functions of the same
order. To meet the convergence condition given in Proposition 6, the parameter in (30) is set to ε = 1.111. The corresponding
feedforward control signals with T = 5, α1 , . . . , α12 , that steer the beam to deform are illustrated in Fig. 3. The evolution of
beam shapes and the regulation error are depicted in Fig. 4. It can be seen that the beam is deformed to the desired shape
and the regulation error vanishes along the entire beam, which confirms the expected performance of the developed control
scheme.
−3

100

x 10

50

−1
−2
−3
−4
0

0

d
w
w d16

α

 d (x), w d (x)
w

0


w d − w d8 1 = 0.0346

w d − w d12 1 = 0.0139

w d − w d16 1 = 0.0037

w d12
w d8
0.2

0.4

−50

α1−12
−100
0.6

x

α1−8

0.8

−150
0

1

α1−16
0.2

0.4

(a)

x

0.6

0.8

1

(b)

Fig. 2. Effect of number of actuators: (a) interpolation accuracy (w
ed (x): desired shape; wdn (x): solution of the steady-state beam equation with n =8, 12
and 16 in-domain actuators); (b) amplitude of steady-state control signals for different setups.

60
α1
α2

40

α3

α1, ..., α12

20

α4
α5

0

α6
α7

−20

α8
α

−40

9

α10
α11

−60

α

12

−80
0

2

4

6

8

10

t

Fig. 3.

Feedforward control signals: α1 , . . . , α12 .

VII. C ONCLUSION
This paper presented a solution to the problem of in-domain control of a deformable beam described by an inhomogeneous
PDE. A relationship between the original model and a system expressed in a standard boundary control form has been

11

(a)
Fig. 4.

(b)

Set-point control: (b) system response; (c) regulation error.

established. Flatness-based motion planning and feedforward control are then employed to explore the degrees-of-freedom
offered by the system, while a closed-loop control is used to stabilize the system around the reference trajectories. The validity
of the developed approach has been assessed through well-posedness and stability analysis. System performance is evaluated by
numerical simulations, which confirm the applicability of the proposed approach. Note that in order to improve the resolution
of manipulation, more actuators may be expected. Nevertheless, for the control of one-dimensional beams, the proposed scheme
contains only one closed loop, allowing a drastic simplification of control system implementation and operation. This is an
important feature for practical applications, such as the control of large-scale deformable micro-mirrors in adaptive optics
systems.
VIII. ACKNOWLEDGMENT
This work is supported in part by a NSERC (Natural Science and Engineering Research Council of Canada) Discovery Grant. The second author is also supported in part by the Fundamental Research Funds for the Central Universities
(#682014CX002EM).
R EFERENCES
[1] Stewart, J. B., Diouf, A., Zhou, Y. Zhou, Bifano, T. G.: ‘Open-loop control of a MEMS deformable mirror for large-amplitude wavefront control’, J. of
the Optical Society of America, 2007, 24, (12), pp. 3827–3833.
[2] Vogel, C. R., Yang, Q.: ‘Modelling, simulation, and open-loop control of a continous facesheet MEMS deformable mirror’, J. of Opt. Soc. Am., 2006,
24, (12), pp. 3827–3833.
[3] Bensoussan, A., Da Prato, G., Delfor, M., Mitter, S. K.: ‘Representation and Control of Infinite-Dimensional Systems’ (Birkhauser, Boston, 2nd ed.,
2006).
[4] Morris, K.: ‘Control of systems governed by partial differential equations’, in Levine, W. S. (Ed): ‘The Control Systems Handbook: Control System
Advanced Methods, Chapter 67’. (CRC Press, 2010, 2nd ed.), pp. 1–37.
[5] Balas, M. J.: ‘Active control of flexible dynamic systems’, Journal of Optimization Theory and Applications, 1978, 25, (3), pp. 415–436.
[6] Guo, B. Z. Wang, J. M. Wang, Yang, K.: ‘Stabilization of an Euler-Bernoulli beam under boundary control and non-collocated observation’, Systems
and Control Letters, 2008, 57, pp. 740–749.
[7] Krstić, M., Smyshyaev, A.: ‘Boundry Control of PDEs: A Course on Backstepping Designs’ (SIAM, New York, 2008).
[8] Luo, Z. H. Luo, Guo, B. Z., Morgul, O.: ‘Stability and Stabilization of Infinite Dimentional Systems with Applications’ (Springer-Verlag, London, 1999).
[9] He, W., Lie, C.: ‘Vibration control of a Timoshenko beam system with input backlash’, IET Control Theory Appl., 2015, 9, (12), pp. 1802–1809.
[10] He, W., Qin2, H., Liu3, J.-K.: ‘Modelling and vibration control for a flexible string system in three-dimensional space’, IET Control Theory Appl., 2015,
9, (16), pp. 2387–2394.
[11] Ammari, K., Tucsnak, M.: ‘Stabilization of Bernoulli-Euler beams by means of a pointwise feedback force’, SIAM J. Control Optimal, 2000, 39, (4),
pp. 1160–1181.
[12] Lasiecka, I., Triggiani, R.: ‘Control Theory for Partial Differential Equations: Continuous and Approximation Theories’ (Cambridge Universioty Press,
Cambridge, 2000).
[13] Ammari, K., Mehrenberger, M.: ‘Study of the nodal feedback stabilization of a string-beams network’, J. Appl Math Comput, 2011, 36, pp. 441–458.
[14] Banks, H. T., Smith, R. C., Wang, Y.: ’Smart Material Structures: Modeling, Estimation and Control’ (John Wiley & Sons, Chichester, 1996).
[15] Chen, G., Delfour, M. C., Krall, A., Payre, G.: ‘Modeling, stabilization and control of serrially connected beams’, SIAM J. Control and Optimization,
1987, 25, (3), pp. 526–546.
[16] Le Gall, P., Prieur, C., Rosier, L.: ‘Output feedback stabilization of a clamped-free beam’, Int. J. of Control, 2007, 80, (8), pp. 1201–1216.
[17] Banks, H. T., del Rosario, R. C. H., Tran, H. T.: ‘Proper orthogonal decomposition-based control of transverse beam vibration: Experimental
implementation’, IEEE Trans. Contr. Syst. Technol., 2012. 10, (5), pp. 717–726.
[18] Meurer, T.: ‘Control of Higher-Dimensional PDEs: Flatness and Backstepping Designs’ (Springer, Berlin, 2013).
[19] Schröck, J., Meurer, T., Kugi A.: ‘Motion planning for piezo-actuated flexible structures: Modeling, design, and experiment’, IEEE Trans. Contr. Syst.
Technol., 2012, 21, (3), pp. 807–819.
[20] Aoustin, Y., Fliess, M., Mounier, H., Rouchon, P., Rudolph, J.: ‘Theory and practice in the motion planning and control of a flexible robot arm using
Mikusinski operators’. Proc. of the Fifth IFAC Symposium on Robot Control, France, Nantes, 1997, pp. 287–293.
[21] Dunbar, W. B., Petit, N., Rouchon, P., Martin P.: ‘Motion planning for a nonlinear Stefan problem’, ESAIM: Control, Optimisation and Calculus of
Variations, 2003, 9, pp. 275–296.
[22] Fliess, M., Martin, P., Petit, N., Rouchon, P.: ‘Active restoration of a signal by precompensation technique’. The 38th IEEE Conf. on Decision and
Control, Phoenix, AZ, Dec. 1999, pp. 1007–1011.
[23] Laroche, B., Martin, P., Rouchon, P.: ‘Motion planning for the heat equation’, Int. J. of Robust Nonlinear Control, 2000, 10, pp. 629–643.
[24] Petit, N., Rouchon, P.: ‘Flatness of heavy chain systems’, SIAM J. Control and Optimization, 2001, 40, (2), pp. 475–495.

12

[25] Petit, N., Rouchon, P., Boueih, J. M., Guerin, F., Pinvidic, P.: ‘Control of an industrial polymerization reactor using flatness’, Int. J. of Control, 2002,
12, (5), pp. 659–665.
[26] Rabbani, T. S., Di Meglio, F., Litrico, X., Bayen, A. M.: ‘Feed-forward control of open channel flow using differential flatness’, IEEE Trans. Contr.
Syst. Technol., 2011, 18, (1), pp. 213–221.
[27] Rudolph, J.: ‘Flatness Based Control of Distributed Parameter Systems’ Shaker-Verlag, Aachen, 2003.
[28] Timoshenko, S., Woinowsky, S.: ‘Theory of Plates and Shells’. (MacGraw-Hill Book Company, New York, 2nd ed., 1959).
[29] Rebarber, R.: ‘Exponential stability of coupled beam with dissipative joints: A frequency domain approach’, SIAM J. Control and Optimization, 1995
33, (1), pp. 1–28.
[30] Pederson, M.: ‘Functional Analysis in Applied Mathematics and Engineering’ (Chapman and Hall CRC, London, 2000).
[31] Adams, R. A., Fournier, J. F.: ‘Sobolev Spaces’ (Academic Press, New York, 2nd ed., 2003).
[32] Abramowitz, M., Stegun, I. A.: ‘Handbook of Mathematical Functions’ (Dover Publications, New York, 1972).
[33] Jacob, B., Zwart, H.: ‘Equivalent conditions for stabilizability of infinite-dimensional systems with admissible control operators’, SIAM journal on control
and optimization, 1999, 37, (5), pp. 1419–1455.
[34] Pazy, A.: ‘Semigroups of Linear Operators and Applications to Partial Differential Equations’ (Applied Math. Sciences 44, Springer, New-York, 1983).
[35] Coron, J. M.: ‘Control and Nonlinearity’ (American Mathematical Society, 2007).
[36] Rodino, L.: ‘Linear Partial Differential Operators in Gevrey Spaces’ (World Scientific, River Edge, 1993).
[37] Yang, B.: ‘Stress, Strain, and Structural Dynamics: an Interactive Handbook of Formulas, Solutions, and MATLAB Toolboxes’ (Elsevier Academic
Press, Burlington, 2005).

IX. A PPENDICES
A. Proof of Proposition 6
We prove the convergence of the power series (28) using Cauchy-Hadamard Theorem. The convergence of (29) then follows
easily using the same argument.
Denote in (27):
bn =

n
X

x4k+1
(4k + 1)!(4(n−k) + 2)!
k=0
!
n
X
x4k+3
(2n)
(−1)n φj (t).
−
(4k + 3)!(4(n−k))!

(38)

k=0

p
n

Then, (28) converges if lim supn→∞ |bn | < 1.
Thus, for 0 < xj < 1 and x ∈ [0, 1], bn can be majorized as:
|bn | <

n
X

1
(4k + 1)!(4(n−k) + 2)!
k=0
!
n
X
1
(2n)
+
φj (t)
(4k + 3)!(4(n−k))!
k=0

n
X

<
+

k=0
n
X

4k + 2
(4k + 2)!(4(n−k))!

k=0

4k + 2
(4k + 2)!(4(n−k))!

≤2(4n + 2)

n
X
k=0

<2(4n + 2)

!
(2n)

φj

(t)

1
(2n)
φj (t)
(4k + 2)!(4(n−k))!

2M 16n ((2n)!)σ
16n
(2n)
φj (t) ≤ 2n
.
(4n + 2)!
K
(4n + 1)!

Therefore,
lim sup

p
n

bn ≤ lim sup 16

n→∞

Applying Stirling’s formula

√
n

1/n

((2n)!)σ/n
((4n + 1)!)1/n
n→∞
2nσ/n
((2n)!)1/(2n)
16
= 2 lim sup
.
K n→∞ ((4n + 1)!)1/(4n+1) (4n+1)/n


2M
K 2n

n! ' (n/e) yields:
lim sup
n→∞

p
n

bn ≤

16e4−2σ
(2n)2nσ/n
lim
sup
.
(4n+1)/n
K2
n→∞ (4n + 1)

(39)

13

Since
(2n)2nσ/n
(4n+1)/n
n→∞ (4n + 1)

σ < 2;

 0,
2σ+ 3σ
n
(2n)
1
= lim sup
, σ = 2;
1 =

n→∞ (4n + 1)4+ n
 16
∞, σ > 2,
lim sup

(40)

we can conclude by Cauchy-Hadamard Theorem that (28) converges for σ < 2 with a radius of convergence of infinity and
for σ = 2 if K 2 > 1. The series (28) diverges if σ > 2.
B. Proof of Claim 7
Consider the first series in (34). Fixing m > 0, for x, xj ∈ (0, 1), j = 1, 2, . . . , N , we have
(−1)k+1 m2k−1 (x − xj )2k+2
m2k−1
≤
.
(2k − 1)
(2k + 2)!
(2k − 1)!
Since limn→∞

1
√
n
n!

= 0, it follows that for any a > 0
r
lim

n→∞

n

1
an
= a lim √
= 0.
n
n→∞
n!
n!

P∞ m2k−1
Thus k=0 (2k−1)!
is convergent. We conclude that the first series in (34) is uniformly convergent. The convergence of the
other terms in (34) can be proved in the same way.
C. Proof of Claim 8
We denoted by [G(xi , ξj )]N ×N the matrix in the left-hand-side of (35) formed by Green’s functions. We argue by contradiction. If, otherwise, [G(xi , ξj )]N ×N is not invertible, then it is of rank less than N . Without loss of generality, assume that
there exist N − 1 constants, k1 , k2 , ..., kN −1 , such that
G(x1 , ξN ) =
G(x2 , ξN ) =

N
−1
X
i=1
N
−1
X

ki G(x1 , ξi ),
ki G(x2 , ξi ),

i=1

..
.
G(xN , ξN ) =

N
−1
X

ki G(xN , ξi ).

i=1

The above equations show that
G(x, ξN ) =

N
−1
X

ki G(x, ξi )

(41)

i=1

has N different positive solutions x1 , x2 , · · · xN , xi ∈ (0, 1), i = 1, . . . , N .
We consider two cases:
(i) If N > 3, since ξi , i = 1, . . . , N, are distinguished, G(x, ξj ), i = 1, . . . , N, are all different from each other. Hence
G(x, ξN ) 6≡ k1 G(x, ξ1 ) + k2 G(x, ξ2 ) + · · · + kN −1 G(x, ξN −1 ).
Note that G(x, ξj ), j = 1, . . . , N, are of order at most 3, then (41) has at most 3 different solutions in R, which is a contradiction.
(ii) If N ≤ 3, it is easy to check that (41) has a solution x = 0, and a pair of solutions x = x0 and x = −x0 near the origin
0. By the assumption, (41) has N different positive solutions, then it must be N = 1, which leads to a contradiction with the
non-invertible property of G(x1 , ξ1 ) 6= 0.

