Planning Trajectories for a Class of Linear Partial Differential Equations: An Introduction
01/10/2017Résumé
Auteurs
Joachim Rudolph |
Métriques
<resource xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://datacite.org/schema/kernel-4" xsi:schemaLocation="http://datacite.org/schema/kernel-4 http://schema.datacite.org/meta/kernel-4/metadata.xsd"> <identifier identifierType="DOI">10.23723/545:2004-1/20067</identifier><creators><creator><creatorName>Joachim Rudolph</creatorName></creator></creators><titles> <title>Planning Trajectories for a Class of Linear Partial Differential Equations: An Introduction</title></titles> <publisher>SEE</publisher> <publicationYear>2017</publicationYear> <resourceType resourceTypeGeneral="Text">Text</resourceType><dates> <date dateType="Created">Sun 1 Oct 2017</date> <date dateType="Updated">Sun 1 Oct 2017</date> <date dateType="Submitted">Mon 18 Feb 2019</date> </dates> <alternateIdentifiers> <alternateIdentifier alternateIdentifierType="bitstream">43535b8eff6267e695b5b654840167dc87713c14</alternateIdentifier> </alternateIdentifiers> <formats> <format>application/pdf</format> </formats> <version>34084</version> <descriptions> <description descriptionType="Abstract"></description> </descriptions> </resource>
1 Planning Trajectories for a Class of Linear Partial Differential Equations: An Introduction Joachim Rudolph Abstract— Trajectory planning for a rather large class of boundary controlled linear distributed parameter systems in one space dimension is considered. Operational calculus and basic algebra leads to operational solutions which are then interpreted by using either power series expansions and Gevrey functions or integral representations. Several examples illustrate the method: models of a tubular chemical reactor, a piezo-electric beam, heat exchangers, and heavy ropes. Index Terms— Boundary control, distributed parameter sys- tems, flatness, motion planning, partial differential equations. I. INTRODUCTION WE consider a class of linear distributed parameter sys- tems in one space dimension. The controls are assumed to be lumped variables, acting at isolated points of the spatial domain. These systems comprise partial differential equations (p.d.e.). They are transformed into ordinary differential equa- tions (o.d.e.) by using the operational calculus introduced by Mikusiński [1], [2]. The solutions of the resulting (ordinary) boundary value problems are operational functions. The cal- culations are the same as with an “innocent” (or “formal”) use of Laplace transform calculus, which supposes the existence of the integrals involved. We are particularly interested in the problem of planning trajectories realizing finite time transitions between two equi- librium profiles (solutions of the stationary boundary value problem). The (open-loop) control is then simply obtained as a consequence of the trajectory planning. We solve this problem by independently assigning trajectories to the components of a “basic output”, also called “flat output” in reminiscence to flat nonlinear finite dimensional control systems [3]–[5]. The trajectories of all other system variables follow from the trajectories of the “flat output”. Such an open-loop control design can readily be combined with stabilizing control laws, which are not treated here. We could employ the algebraic framework of module the- ory. The methods presented here are then generalizations of those derived when dealing with so-called π-free linear delay systems [6], [7] (see also [8] for the concept of π- flatness of nonlinear delay systems). The “basic outputs”, or “flat outputs”, are then defined as bases of modules. We are not discussing this mathematical framework here, but references like [7], [9], and in particular [10], give all the details about this approach in a manner which is consistent with the approach followed here. Manuscript submitted January 31, 2003. J. Rudolph is with the Institut für Regelungs- und Steuerungstheorie, TU Dresden, 01069 Dresden, Germany The present introduction to the subject is structured as follows. In the next section we introduce the class of sys- tems considered. In Section III we discuss their operational solutions in order to introduce the notion of basic (or flat) outputs in Section IV. Using the parameterization of the trajectories with such a basic output allows us to reduce the trajectory planning problem to the time domain interpretation of operators and operational functions. This is the subject of Section V. Sections VI to IX serve to illustrate the method on examples: a chemical reactor modeled as a parabolic equation, a piezo-electric beam model, and models of heat exchangers and heavy ropes. A short discussion of relevant references (necessarily incomplete) is given in conclusion — see also [11], [12]. II. CLASS OF SYSTEMS CONSIDERED Some of the variables of the systems under consideration represent “distributed quantities”, w1, . . . , wp, which depend both on time t and a single spatial variable x, the others represent “lumped quantities”, v1, . . . , vr, depending on time only. We are, thus, given homogeneous systems of linear p.d.e. for w = (w1, . . . , wp): α X j=0 β X k=0 p X l=1 ai,j,k,l(x) ∂j+k ∂xj∂tk wl(x, t) = 0, (1) x ∈ Ωi = [xi,0, xi,1] ⊂ R, t ≥ 0, i = 1, . . . , q, with (appropriate) boundary conditions α X j=0 β X k=0 p X l=1 bi,j,k,l ∂j+k ∂xj∂tk wl(x, t) + β̄ X k̄=0 r X l̄=1 ci,k̄,l̄ dk̄ dtk̄ vl̄(t) = 0, (2) x = xi ∈ Γ = {x1,0, . . . , xq,0, x1,1, . . . , xq,1}, t ≥ 0, i = 1, . . . , α̃, as well as (appropriate) initial conditions ∂k wl ∂tk (x, 0) = w0,l,k(x), x ∈ Ωl ∈ {Ω1, . . . , Ωq}, 0 ≤ k ≤ γl, l = 1, . . . , p, dk vl dtk (0) = v0,l, l = 1, . . . , r, 0 ≤ k ≤ µl. (3) The domains of definition of the variables wl, l = 1, . . . , p, are closed intervals Ωl ⊂ R. Thus, in (1), the coefficients ai,j,k,l(x) are different from zero only if the wl are defined 2 on Ωi, i.e., if Ωl = Ωi. The same holds true for the coefficients bi,j,k,l in the boundary conditions (2): Condition i concerns a boundary point xi of one or more intervals Ωl, and bi,j,k,l is different from zero only if xi is a boundary point of Ωl. If the coefficients ai,j,k,l in (1) depend on x we suppose they are real valued functions which are analytic on the whole interior of an interval, else they are real numbers. Coefficients bi,j,k,l and ci,k̄,l̄ in the boundary conditions (2) are taken from R, too. Notice that we consider only homogeneous p.d.e. They involve variables defined on (one or more) finite intervals, which can share common boundary points. The boundary conditions (2) may comprise o.d.e. in (lumped) variables vi, i = 1, . . . , r. In this case, the coefficients bi,j,k,l are identically zero for some i ∈ {1, . . . , α̃}, and j, k, l. Control variables are chosen between the lumped variables vi, i = 1, . . . , r, which means we do not treat distributed control. A (lumped) control acting at an interior point x1 of an interval [x0, x2] is considered by introducing two intervals [x0, x1] and [x1, x2]. Doing so, the control becomes a bound- ary control and variables are defined on both intervals. Finally, we consider only solutions satisfying the p.d.e. on the whole domain of definition, including the boundary. The ordinary boundary value problem α X j=0 p X l=1 ai,j,0,l(x) dj dxj w̃l(x) = 0, x ∈ Ωi = [xi,0, xi,1] ⊂ R, i = 1, . . . , q, α X j=0 p X l=1 bi,j,0,l dj dxj w̃l(x) + r X l̄=1 ci,0,l̄ ṽl̄ = 0, x = xi ∈ Γ, t ≥ 0, i = 1, . . . , α̃, resulting from (1) and (2) for variables which do not depend on time are called stationary profiles. They depend on stationary values ṽj, j = 1, . . . , r, of the lumped variables vj, j = 1, . . . , r. If (w1, . . . , wp, v1, . . . , vr) is a solution of the prob- lem defined by (1), (2), and (3), and if (w̃1, . . . , w̃p, ṽ1, . . . , ṽr) is a solution of the associated ordinary boundary value prob- lem, then — as a consequence of the linearity of the equations — (w1 +w̃1, . . . , wp +w̃p, v1 +ṽ1, . . . , vr +ṽr) is a solution of (1), (2), too. In other words, choosing zero initial conditions does not mean any restriction for the trajectory planning problem we consider, which consists in a transition between two stationary profiles. III. SOLUTION VIA OPERATIONAL CALCULUS We now replace derivation w.r.t. t by the operator s (in the sense of Mikusiński’s operational calculus1 ), using the fact that the initial conditions have been taken to be zero. (More exactly, one has ˆ ϕ̇ = sϕ̂ − ϕ(0) with ϕ(0) = 0.) This way we 1With the Laplace transform [13] we analogously get ŵ(x, s) = R ∞ 0 e−stw(x, t)dt and û(s) = R ∞ 0 e−stu(t)dt. get the o.d.e. system α X j=0 β X k=0 p X l=1 ai,j,k,l(x)sk dj dxj ŵl(x) = 0, (4) x ∈ Ωi = [xi,0, xi,1] ⊂ R, i = 1, . . . , q, with boundary conditions α X j=0 β X k=0 p X l=1 bi,j,k,lsk dj dxj ŵl(x) + β̄ X k̄=0 r X l̄=1 ci,k̄,l̄sk̄ v̂l̄ = 0, (5) x = xi ∈ Γ = {x1,0, . . . , xq,0, x1,1, . . . , xq,1}, i = 1, . . . , α̃, where s is a parameter. Here ŵl : Ωl → M, l = 1, . . . , p, denotes the operational function associated with wl, the image of which for every fixed x ∈ Ωl is an operator, and v̂j, j = 1, . . . , r, are the operators associated with the functions vj. The definition of the functions ŵl, l = 1, . . . , p, can be extended on Ω = ∪l=1,...,pΩl by posing ŵl(x) = 0 on Ω\Ωl. Thus, ŵ can be interpreted as a vector valued function defined on Ω and the system of o.d.e. (4) can be written as α X j=0 Ai,j(x) dj dxj ŵ(x) = 0, x ∈ Ω, i = 1, . . . , q, (6) the boundary conditions (5) as α X j=0 Bi,j dj dxj ŵ(x) + Ci v̂ = 0, x = xi ∈ Γ, i = 1, . . . , α̃. (7) Here Ai,j(x), i = 1, . . . , q, j = 0, . . . , α, are matrices (of dimension 1 × l) over the ring of polynomials in s with coefficients in the ring of real analytic functions on appropriate intervals Ωi, and Bi,j, i = 1, . . . , α̃, j = 0, . . . , α, and Ci, i = 1, . . . , α̃, are matrices (with dimensions 1 × l and 1 × r, respectively) over the ring of polynomials in s with coefficients in C. Remark: Notice that the relations between the components of v could as well represent a time invariant linear system with (constant) delays or that delayed variables could occur in the boundary conditions. The components of the matrices Ci, i = 1, . . . , α̃, would then simply be operators in the ring C[s, e−τ1s , . . . , e−τρs ]. Assume we know a system of solutions of system (6) which allows us to satisfy the boundary conditions, i.e., an appropriate set of M-linearly independent solutions Ŵi = (ŵi,1, . . . , ŵi,p), i = 1, . . . , n, with components ŵi,l, i = 1, . . . , n, l = 1, . . . , p, which are operational functions on Ω. Then their linear combinations ŵ(x) = n X i=1 Ki Ŵi(x), Ki ∈ M, x ∈ Ω (8) are solutions of (6) (with values in M), too [2]. In case the coefficients ai,j,k,l in (4) are constants, the well-known classical methods allow us to calculate a fundamental system 3 of solutions such that (8) provides all solutions. One calculates the roots λi ∈ M of the characteristic equation and composes the solution with components of the type xκi−j eλix , j = 1, . . . , κi, according to the multiplicity κi of the λi. Here, only such roots λi must be considered which are logarithmic in the sense of Mikusiński, which means that the associated operational functions eλix exist (for instance eisx , where i = √ −1, does not exist [2]). If the coefficients in the o.d.e. depend on the spatial variable x, then obtaining the solution is much more difficult in general. Sometimes, however, one may use transformations of the variables (see the examples of the heavy ropes in Section IX and of the flexible beam in [14]) or power series in x [15]. Remark: For p.d.e. with constant coefficients, dealing with roots ±λ it is often useful to use hyperbolic or trigonometric rather than exponential functions. This simplifies the satisfaction of the boundary conditions. Moreover, we will use operational functions which admit direct interpretations as functions of time, with forward or backward shift (or “delay”) operators, eαs , α ∈ R, or power series in s, for instance. The examples discussed in the sequel will illustrate these techniques. Using formula (8) for ŵ in the boundary conditions (7), with (8), leads to a non-homogeneous system of linear equations for the coefficients Kk, k = 1, . . . , n: α X j=0 n X k=1 Bi,j dj Ŵk dxj (xi) Kk + Ci v̂ = 0, i = 1, . . . , α̃, n X k=1 Ŵk(x) Kk = ŵ(x). Arranging the equations appropriately, this system can be rewritten in matrix form: 0 S W K = C v̂ D v̂ ŵ , (9) where S is an (n × n) matrix with entries Si,k = α X j=0 Bi,j dj Ŵk dxj (xi), for k = 1, . . . , n and a choice of n indices from {1, . . . , α̃}, C is a ((α̃−n)×r) matrix, D one (n×r), the rows of which correspond to Ci, i = 1, . . . , α̃, and W = (Ŵ1(x), . . . , Ŵn(x)). Hence, S is a matrix the entries of which are operators (from the field M), while C and D are matrices over C[s]. In contrast, the components of W are operational functions. The equation D v̂ = 0 describes the relations between the lumped variables, which have been combined for further discussion. (If delays occur, as mentioned in the remark on page 2, C and D are defined over C[s, e−τ1s , . . . , e−τρs ].) If the boundary conditions have been chosen appropriately (which will be the case in practical problems), and if the solution formula (8) is rich enough, the ordinary boundary value problem (6), (7) can be solved. Then the matrix S is invertible (α̃ = Pp l=1 nl), and the coefficients in (8) are given by the formula K = S−1 D v̂. Hence, the solution of our problem has the form ŵ = WS−1 D v̂, (10) with v̂ satisfying the equation C v̂ = 0 (often C = 0). The equivalent equations C v̂ = 0 (det S) ŵ − W(adj S)D v̂ = 0 can be written with coefficients from the ring generated by the entries of the matrices W, S, C, and D over C (here adj S is the adjoint of S, which satisfies S−1 = (adj S)/(det S)). We, thus, have linear equations for the system variables, or (more exactly) for the family ŵ = (ŵ1, . . . , ŵp) of operational functions and the family v̂ = (v̂1, . . . , v̂r) of operators, which are involved in the boundary conditions. These equations are at the base of our subsequent discussions. IV. SYSTEMS AND “BASIC OUTPUTS” In the preceding section we have obtained a system of linear equations in a matrix form, linking ŵ1, . . . , ŵp and v̂1, . . . , v̂r: C v̂ = 0 (11a) (det S) ŵ − W(adj S)D v̂ = 0, (11b) where the matrices S, C, and D are composed of operators (from M), C and D even polynomials from C[s], while the entries of matrix W are operational functions. Let S and W be the families of the entries of S and W, and R = C[s, S, W] (the generators of which are not necessarily independent). We, thus, have R-linear relations between the variables, which allows us to identify the system with an R- module [7], [9], [10] — we will not discuss these aspects here. For the moment, let C = 0 in (11a), the family v̂ forms an independent input then. Observe that using another family of variables, b̂ = (b̂1, . . . , b̂r), we can write ŵ = W(adj S)D b̂ (12a) v̂ = (det S) b̂ (12b) without violating the system equation (11b); we have (det S) ŵ − W(adj S)D v̂ = 0. Note that b̂, as the input v̂, represents lumped variables. In analogy to the notion of π-freeness of linear delay systems, see [6], [7] for instance, and with differential flatness of nonlinear finite dimensional systems [3]–[5], the family b̂ is called π-basic, or (det S)-basic, or also π- or (det S)-flat output of the system. For simplification, we also just say flat or basic output, even if det S 6∈ C. In the module theoretic approach to the class of systems under consideration, proposed in [7], [9] and detailed in [10], b̂ is a basis of a module associated with the system and with π = det S. 4 If C 6= 0 in (11a) this equation introduces relations between components of v̂. However, choosing an appropriate matrix N allows us to get ŵ = W(adj S)DN b̂ (13a) v̂ = (det S)N b̂ (13b) where the entries of N belong to C[s, S]. (First, one constructs a π̃-basic (π̃-free) output b̃ of the C[s, S]-module generated by v̂, with v = π̃−1 Mb̃, whence N = π̃−1 M. Accordingly, the system is then called (π̃ det S)-free (or -flat), and b̂ is called (π̃ det S)-basic (or -flat) output — see [10] for further details.) The entries of N belonging to C[s, S] the coefficients in (13) are composed of operators and operational functions from R. Therefore, in order to give a time domain interpretation of the parameterization by a π-basic output it is enough to know the interpretation of these operators, which are introduced when solving the ordinary boundary value problem. It is not required to interprete π−1 ! Remark: The systems considered are not necessarily controllable. If the matrices P = W(adj S)D and (det S)I have a (non-unimodular) common left divisor (over R) there is a non-controllable (namely an autonomous) part. However, this non-controllable part does not play any role in the problem of planning trajectories starting in a stationary profile, because this part will “rest at zero”. Example: An example of trajectory planning for a non- controllable system is studied in Section IX-C, where we study the planar motion of two heavy ropes of equal length the suspension point of which can be displaced horizontally. V. TRAJECTORY PLANNING AND CALCULATION OF CONTROLS We will study two cases in the sequel, on the one hand oper- ators and operational functions admitting an interpretation via series expansions depending on derivatives of the flat outputs, on the other hand operators and operational functions leading to compact support functions, which leads to distributed finite delays and advances. Let P : Ω → M be an operational function admitting a power series expansion P n≥0 pn(x)sn , x ∈ Ω (see [1], [2]). Then, the formula ŵ = P b̂ leads to a series P n≥0 pn(x)sn b̂, and, if the derivatives b(n) are all zero at t = 0 and for n > 0, then we have w(x, t) = X n≥0 pn(x) b(n) (t). Convergence of this series will be discussed case by case here. To this end, we will use functions (of t) all derivatives of which are zero at the beginning, t = 0, and at the end, t∗, of a time interval [0, t∗] used for the transition. These functions will be Gevrey functions of a class (or order) α > 1, with α ≤ 2 in our examples, which means sup τ∈R+ |b(n) (τ)| ≤ m (n!)α γn , n ≥ 0, where α, m, and γ are constant real parameters [16], [17]. Remark: 1. Notice that if functions b are not identically zero they are not analytic at t = 0, since an analytic function all derivatives of which are zero at a point are constant. Gevrey functions of class α = 1 are analytic, whence α > 1 is required here. 2. The vector valued function P and the operators in S are treated similarly. An example of an appropriate Gevrey function is ηd : R → R with ηd(t) = 0, for t < 0, 1 2 + 1 2 tanh 2(2 t t∗ − 1) (4 t t∗ (1 − t t∗ ))d ! , for t ∈ [0, t∗], 1, for t > t∗ (14) where d ≥ 1 is a real parameter. It is a Gevrey function of order (1 + d)/d. It is smooth and its derivatives satisfy η̃ (n) d (0) = η̃ (n) d (1) = 0, n > 0. (Successive derivatives can be calculated by recursion [10] so that numerical implementation is easy.) The choice of the parameter d allows us to vary the trajectories used for the transition (see Fig. 1). 0 1 0 1 y t d=1.1 d=10 d = 10 d = 1.1 t/t∗ ηd Fig. 1. Graph of the function ηd in (14) for two different values of the parameter d. Example: The use of power series is illustrated on the examples of Sections VI and VII. Now let P : Ω → M be an operational function such that for all x ∈ Ω the corresponding time function p = {P(x, t)} has compact support in I = [t0, t1] ⊂ R+, i.e., p(x, t) = 0 for t 6∈ I, x ∈ Ω. (We then speak about compact support operators2 .) Then interpretation of formula ŵ = P b̂ as a convolution product yields w(x, t) = Z t 0 p(x, τ) b(t − τ) dτ = Z t1 t0 p(x, τ) b(t − τ) dτ. We see that for a given time instant T the values of b on the interval [T −t1, T −t0] are required when calculating w(x, T), 2These operators correspond to those which with Laplace calculus satisfy the conditions of a Paley-Wiener Theorem. 5 in other words, values from the past. Hence, we have a (finite) distributed delay. Often also shift operators ehs with h ∈ R+ occur. They give rise to distributed advances (predictions). Remark: Notice that, in contrast to operators leading to series, in the compact support case functions defining trajectories of b need not be smooth. Example: The compact support operators occur in the heat exchanger examples of Section VIII and the heavy rope example in Section IX. VI. EXAMPLE: AN ISOTHERMAL TUBULAR REACTOR We consider the linear parabolic equation with constant coefficients λ0 ∂2 w(x, t) ∂x2 − ν0ρ0 ∂w(x, t) ∂x + µ0w(x, t) = ρ0 ∂w(x, t) ∂t , −L ≤ x ≤ 0, t ≥ 0. (15) This p.d.e. can serve as a model of a chemical or biotechno- logical reactor in which diffusion effects are not negligible in comparison to the forced convection, the latter being supposed of plug flow type. In such a linearized model µ0w can represent the reaction effects. Then, concentration of a species w depends on time and on the axial coordinate of the reactor x. Dependency in the other directions is neglected. The constant parameters are λ0 > 0 for diffusion, ν0 ≥ 0 for convection, ρ0 > 0 for density, and µ0 for reaction. The model includes the heat equation as a special case (with ν0 = µ0 = 0). Generalizations to varying coefficients, to systems of equations, and beyond the linear case can be found in [15], [18]–[22]. The boundary conditions can be written ∂w ∂x (0, t) = 0, t ≥ 0 (16a) w(−L, t) = u(t), t ≥ 0. (16b) Here u(t) plays the role of the control variable; the more realistic case ν0ρ0w(−L, t) − λ0 ∂w ∂x (−L, t) = v(t), t ≥ 0, can be treated similarly. As soon as a series representation has been found for w the input v will be obtained by evaluating the series at x = −L. The boundary condition involving the input does not play any role in the representation of w for this parabolic equation, as will be seen in the sequel. Being interested in transitions between stationary profiles, initial conditions can be taken to be zero: w(x, 0) = 0, t ≥ 0. A. Introducing a basic output Introducing the operator s for the differentiation w.r.t. t (with zero initial conditions) we get the o.d.e. λ0 d2 ŵ dx2 (x) − ν0ρ0 dŵ dx (x) + µ0 ŵ(x) − ρ0 s ŵ(x) = 0, −L ≤ x ≤ 0, (17) with the boundary conditions dŵ dx (0) = 0 and ŵ(−L) = û. The characteristic equation λ0 p2 − ν0ρ0 p + (µ0 − ρ0s) = 0 has roots p1 = β0 + β and p2 = β0 − β, where β0 = − ν0ρ0 2λ0 and β = q ν2 0 ρ2 0 − 4λ0(µ0 − ρ0s). In view of the boundary conditions we determine a solution of the form ŵ(x) = eβ0x K1 cosh βx + K2 sinh(βx) β . Taking eβ0x (sinh(βx))/β as one of the solutions of the fundamental system, instead of eβ0x sinh(βx), guarantees us that the two operational functions, which will occur in the coefficients, admit power series expansions in s, since β2 has the form as + b. The boundary condition at x = 0 yields K2 = −β0 K1, and at x = −L we get û = e−β0L cosh(βL) + β0 sinh(βL) β K1, which allows us to determine the parameter K1. For the solution it follows ŵ(x) = P(x) P(−L) û with P(x) = eβ0x cosh(βx) − β0 sinh(βx) β . (18) The system, thus, reads Q ŵ − P û = 0 with Q = P(−L). Therefore, choosing ŵ = P b̂ and û = Q b̂ equation Q ŵ − P û = 0 will be fulfilled. We have obtained a parameterization of the solution with a Q-basic output b̂. Evaluating the operational function ŵ at x = 0 yields b̂ = ŵ(0). The variable b corresponds to the basic output which thus admits an interpretation as the concentration at the reactor outflow. 6 B. Planning trajectories Developing the hyperbolic functions in the operational func- tion P in (18) in their power series in the argument (βx) (which depends on s) yields ŵ(x) = eβ0x X n≥0 (1 − β0) (βx) 2n (2n)! b̂ = e− ν0ρ0 2λ0 x × X n≥0 1 + ν0ρ0 2λ0 (x) 2n (ν2 0 ρ2 0 − 4λ0(µ0 − ρ0s))n (2n)! b̂, i.e., a power series in s. This series can be interpreted by replacing s by differentiations w.r.t. time t (assuming zero initial conditions): w(x, t) = e− ν0ρ0 2λ0 x × X n≥0 1 + ν0ρ0 2λ0 (x) 2n (ν2 0 ρ2 0 − 4λ0(µ0 − ρ0 d dt ))n (2n)! b(t). (19) The series for the control u(t) follows from the latter by evaluating it at x = 0. Choice of the trajectory: Let b : R → R be a C∞ -function, with b(t) = 0 for t < 0 and the following properties: • b(n) (0) = 0, n ≥ 0, for t = 0. • The derivatives are bounded as follows: sup τ∈R+ |b(n) (τ)| ≤ m (n!)α γn , with α < 2, n ≥ 0, and m, γ ∈ R, which means b is a Gevrey function of class α < 2. Then, the series in eq. (19) converges, and the pair (w(x, t), u(t)) forms a solution of the trajectory planning problem considered. Sketch of a convergence proof [20] (compare with [18], [23]): For discussing convergence it is useful to use another representation of the series, namely a power series in x: w(x, t) = X k≥0 ak(t) xk k! . Substituting this formula into the p.d.e. (15) leads to a recur- sion formula for the coefficients ak, k ≥ 2: ak = 1 λ0 (ρ0ȧk−2 + ν0ρ0ak−1 − µ0ak−2). (20) The boundary condition at x = 0 gives a1 = 0, and evaluating (20) at x = 0 yields a0(t) = w(0, t) = b(t). By an induction one can show that for l, k ≥ 0 one has sup t∈R+ |a (l) k (t)| ≤ mMk γl ((l + k)!)α (k!)α/2 (21) with constants M, γ, and m and α < 2. For k = 0 this follows from choosing a Gevrey function of class α < 2 for y = a0. Let (21) hold true for 1, . . . , k − 1. Differentiating (20) for k ≥ 2 it follows |a (l) k | ≤ 1 λ0 (ρ0|a (l+1) k−2 | + ν0ρ0|a (l) k−1| + |µ0||a (l) k−2|), whence sup t∈R+ |a (l) k | ≤ ρ0mMk−2 ((l + k − 1)!)α λ0γl+1((k − 2)!)α/2 + ν0ρ0mMk−1 ((l + k − 1)!)α λ0γl((k − 1)!)α/2 + |µ0|mMk−2 ((l + k − 2)!)α λ0γl((k − 2)!)α/2 ≤ mMk ((l + k)!)α γl(k!)α/2 1 λ0 ρ0(k(k − 1))α/2 γM2(l + k)α + ν0ρ0kα/2 M(l + k)α + |µ0|(k(k − 1))α/2 ((l + k)(l + k − 1))αM2 ≤ mMk ((l + k)!)α γl(k!)α/2 1 λ0 ρ0 γM2 + ν0ρ0 M + |µ0| M2 . The expression in parentheses in the line before last is com- posed of a bounded expression and two which are monotoni- cally decreasing with l and k. From this, we get the last line bracket. Taking this bracket equal to λ0 leads to a quadratic in M, the larger root of which gives M. This means (21) holds true. The radius of convergence of the series follows from the Cauchy-Hadamard formula3 . With (21), for l = 0 and α < 2, we get R = limk→∞ inf t∈R+ k! |ak(t)| 1 k = lim k→∞ (k!)(2−α)/(2k) M = ∞. Remark: A (classical) transformation of the dependent variable w transforms equation (15) into the heat equation — see also [23]. Proving convergence is then simple by the Cauchy-Hadamard formula [4]. The direct approach used here is interesting because it can be generalized beyond the linear constant coefficient case [20]–[22], [24]. Systems of equations of this type, modeling balances for several species in the reactor, are treated in [19]. VII. EXAMPLE: A PIEZO-ELECTRIC BEAM We consider the model of a piezo-electric beam, as sketched in Fig. 2. The beam is composed of two connected laminated beams of piezo-electric material, with thickness d, n layers each. The electrodes cover the whole surface. The beam, clamped at one end, has length l and width w. Applying a voltage −V + ≤ V ≤ V + leads to a planar bending (x1, x3). 3Inspecting the proof one sees that α = 2 leads to a finite radius of convergence R ≥ 1/M. 7 x 1 x 3 u 3 u 1 Fig. 2. Sketch of a piezo-electric beam The (small) deformations u3 in direction of x3 are reasonably modeled by the Euler-Bernoulli beam equation: µ ∂2 u3 ∂t2 (x1, t) + B ∂4 u3 ∂x4 1 (x1, t) = 0, 0 ≤ x1 ≤ l, t ≥ 0, (22) with boundary conditions (for t ≥ 0): u3 (0, t) = 0, ∂u3 ∂x1 (0, t) = 0, B ∂2 u3 ∂x2 1 (l, t) + kV (t) = 0, ∂3 u3 ∂x3 1 (l, t) = 0. (23) The constant parameters k, µ, and B are given in [25]. A. Introduction of a basic output It is useful to normalize with a new spatial variable x = x1/l and a normalized time τ = √ B l2√ µ t, and to introduce w(x, τ) = u3(x1, t) for the deformation. Then, with ∂ ∂x1 = 1 l ∂ ∂x and ∂ ∂t = √ B √ µl2 ∂ ∂τ one gets the dimensionless model ∂4 w ∂x4 (x, τ) + ∂2 w ∂τ2 (x, τ) = 0, x ∈ [0, 1], τ ≥ 0, with boundary conditions w(0, τ) = 0, ∂w ∂x (0, τ) = 0, ∂2 w ∂x2 (1, τ) = −u(τ), ∂3 w ∂x3 (1, τ) = 0 involving the dimensionless control u(τ) = kl2 B V (τ). Initial conditions are w(x, 0) = 0, ∂w ∂τ (x, 0) = 0. With the operator s we get the o.d.e. s2 ŵ(x) + d4 ŵ dx4 (x) = 0, x ∈ [0, 1] with boundary conditions ŵ (0) = 0, dŵ dx (0) = 0, d2 ŵ dx2 (1) = −û, d3 ŵ dx3 (1) = 0. The solution can be given in the form ŵ(x) = A cosh √ isx + B sinh √ isx + C cos √ isx + D sin √ isx (where i = √ −1), and the parameters A, B, C, and D satisfy the linear system of equations A + C = 0 B + D = 0 A is cosh √ is + B is sinh √ is − C is cos √ is − D is sin √ is = −û A sinh √ is + B cosh √ is + C sin √ is − D cos √ is = 0 corresponding to the boundary conditions. Eliminating the parameters A and B leads to sinh √ is − sin √ is C + cosh √ is + cos √ is D = 0 is cosh √ is + cos √ is C + is sinh √ is + sin √ is D = û cos √ isx − cosh √ isx C + sin √ isx − sinh √ isx D = ŵ(x). The relations between û, ŵ, C, and D read S1,1 S1,2 S2,1 S2,2 P1 P2 C D = 0 û ŵ . The equation S1,1C + S1,2D = 0 is satisfied if C = −S1,2b̂, D = S1,1b̂, whence (S1,1S2,2 − S1,2S2,1) b̂ = Q b̂ = û (P2S1,1 − P1S1,2) b̂ = P b̂ = ŵ and ŵ = PQ−1 û. (An appropriate choice of the boundary conditions implies Q 6= 0, or more generally the existence of Q−1 .) We proceed in a similar fashion if ŵ and û comprise several components. We are, thus, given a system representation Q ŵ(x) = P(x) û, with P(x) = cosh √ isx − cos √ isx cosh √ is + cos √ is − sinh √ isx − sin √ isx sinh √ is − sin √ is (24a) Q = −is 2 + cosh √ is cos √ is . (24b) Moreover, this calculation implies that ŵ = P b̂, û = Q b̂, 8 i.e., b̂ may be used as a Q-basic (or Q-flat) output. With the “addition theorems” for the hyperbolic functions we can rewrite P and Q in such a way that the formulae in (24) read û = − h 2 + cosh √ is(1 + i) + cosh √ is(1 − i) i b̂ (25) and ŵ = h cosh √ isi (1 − x) − cosh √ is (1 − x) + 1 + i 2 cosh √ is (1 + ix) − cosh √ is (i − x) + 1 − i 2 cosh √ is (1 − ix) − cosh √ is (i + x) i i s b̂. Developing the hyperbolic cosine in its power series allows us to conclude for ŵ(x) = 2 ∞ X k=0 (−1) k (4k + 2)! h (1 − x) 4k+2 + Re n (x + i) 4k+2 o +Im n (x − i) 4k+2 oi s2k b̂, and the corresponding formula for û follows either by differ- entiating twice w.r.t. x and evaluating the result or directly from formula (25): û = −2 " 1 + ∞ X k=0 1 (4k) ! (2s)2k # b̂. B. Planning trajectories The series of the deflection trajectory results from the operational representation by interpreting the operator s as time derivation (with zero initial conditions): w(x, t) = 2 ∞ X k=0 (−1) k (4k + 2)! h (1 − x) 4k+2 + Re n (x + i) 4k+2 o + Im n (x − i) 4k+2 oi b(2k) (t), (26) and for the control we have u(t) = −2 " b(t) + ∞ X k=0 4k (4k)! b(2k) (t) # . As in the preceding section it suffices now to choose a trajectory with Gevrey class less than 2 all derivatives of which vanish at t = 0 (b(n) (0) = 0, n ≥ 0) in order to get a converging series and a solution (w(x, t), u(t)) of the problem. Sketch of a convergence proof: The terms ak, k ≥ 0, of the series (26) satisfy |ak| ≤ 6b(2k) zk+1 /(4k + 2)! =: 6ckz, with z = (x + 1)4 . With The Cauchy-Hadamard Theorem and Stirling’s Formula one obtains r = 1 limk→∞ k p | ck | ≥ lim k→∞ γ2k ((4k + 2)!) m((2k)!)α 1/k = lim k→∞ 16γ2 4−α e2α−4 k4−2α = ∞ for the radius of convergence r of the series P k≥0 ckzk if α < 2. Remark: Further details on piezo-electric beams can be found in [26]. In this article the operators for the Euler-Bernoulli type model are deduced from those of the more general Timoshenko beam model. Loads at the free end are treated, as well as beams with several piezo actuators each covering a part of the beam. In contrast to the Euler-Bernoulli beam models Timoshenko beam models do not lead to power series but to (finite) distributed delays and advances [26], [27], similar to those discussed in Sections VIII and IX. t [ s ] t i p d e f l e c t i o n 1 0 - 4 [ m ] 1 . 2 1 . 0 0 . 8 0 . 6 0 . 4 0 . 2 0 . 0 - 0 . 2 0 . 0 - 0 . 0 3 0 . 0 1 - 0 . 0 2 - 0 . 0 1 0 . 0 2 0 . 0 3 0 . 0 4 0 . 0 5 0 . 0 6 c a l c u l a t i o n m e a s u r e m e n t measured calculated m ) / (s) t / ( w Fig. 3. Experimental result: bending of a piezo-electric beam. Figure 3 illustrates the usefulness of the approach. We see here the result of an experiment made in cooperation with W. Haas at Johannes Kepler University of Linz (Austria) [25]. The control computed using the series has been used as a reference signal of a feedback loop. Another experimental result based on a beam model and the series approach, for a flexible robot arm, can be found in [14]. VIII. EXAMPLE: HEAT EXCHANGERS The simplest type of heat exchangers consists in two (cylin- dric) tubes mounted one into the other (with common axis of symmetry), cf. Fig. 4. x = 0 Ti(x, t) ve Te(x, t) vi x = 1 Ti(0, t) u(t) = Te(0, t) Fig. 4. Sketch of a heat exchanger. For the sake of simplicity, assume the thermal exchange between the outer tube and the environment being negligible, and neglect also heat conduction. The mathematical model of 9 such a typical heat exchanger can then be written (cf. [28]): ∂Ti ∂t + vi ∂Ti ∂x = a (Te − Ti) (27a) ∂Te ∂t + ve ∂Te ∂x = b (Ti − Te) . (27b) In these equations Te(x, t) and Ti(x, t) denote the temperature fields in the exterior and the interior tube, respectively, with x ∈ [0, 1] the spatial variable, and time t ≥ 0. Constant real parameters of the model are: velocities vi and ve, densities ρi, ρe, and specific heat capacities ci, ce of the inner and outer fluid, heat exchange coefficient αi, and radii ri and re. With these physical parameters we get a = 2 αi ci ρi ri and b = 2 αi ri ce ρe (r2 e − r2 i ) . For simplifying notations we consider a unit (normalized) length exchanger. Initial conditions are taken zero: Ti(x, 0) = Te(x, 0) = 0. In both cases we get the same p.d.e., if we use negative interior velocity, vi < 0, for the counter-flow case. However, the fluid enters at x = 1 in this latter case, so that the boundary conditions differ from the parallel to the counter-flow case. Considering the inflow temperature of the exterior tube (at x = 0) as the control, u(t) = Te(0, t), the boundary conditions read Ti(0, t) = Ti0(t) = 0, Te(0, t) = u(t) for the parallel flow, and Ti(1, t) = Ti1(t) = 0, Te(0, t) = u(t) for the counter-flow case. It is useful to divide equations (27) by the velocities and to consider mean and difference parameters: ν = 1 2 1 vi + 1 ve , 4ν = 1 2 1 vi − 1 ve , β = 1 2 a vi + b ve , 4β = 1 2 a vi − b ve . Then equations (27) can be written ∂Ti ∂x + (ν + 4ν) ∂Ti ∂t − (β + 4β) (Te − Ti) = 0 ∂Te ∂x + (ν − 4ν) ∂Te ∂t − (β − 4β) (Ti − Te) = 0. Differential equations of the counter-flow exchanger can be written with positive vi by utilizing those of the parallel flow case and exchanging β and −4β as well as ν and −4ν. A. Operational model and basic output Introducing the operator s for differentiation w.r.t. time and operational functions T̂i(z, s) and T̂e(z, s) leads to the o.d.e. system dT̂i dx + (ν + 4ν) s T̂i − (β + 4β) T̂e − T̂i = 0 dT̂e dx + (ν − 4ν) s T̂e − (β − 4β) T̂i − T̂e = 0. Its characteristic equation (in p) reads p + ν s + β 2 − (4ν s + 4β) 2 − β 2 + 4β2 = 0, with roots p1 = −(ν s + β) + q 4ν2 s2 + 2 4ν 4β s + β 2 = −γ0 + γ1 p2 = −(ν s + β) − q 4ν2 s2 + 2 4ν 4β s + β 2 = −γ0 − γ1 where γ0 = (ν s+β) and γ1 = q 4ν2 s2 + 2 4ν 4β s + β 2 . The boundary conditions are T̂i(0) = 0, T̂e(0) = û for the parallel flow, T̂i(1) = 0, T̂e(0) = û for the counter-flow case. In view of the boundary conditions we consider a solution of the form T̂i(x) T̂e(x) ! = K1 sinh(γ1x) γ1 + K2 cosh(γ1x) e−γ0x K3 sinh(γ1x) γ1 + K4 cosh(γ1x) e−γ0x . Taking the p.d.e. into account at the boundaries we get an equation of the type SK = Dû for the coefficients K1, . . . , K4, in the parallel flow case: 0 1 0 0 1 (s 4ν + 4β) 0 −(β + 4β) 0 0 0 1 0 −(β + 4β) 1 −(s 4ν + 4β) K1 K2 K3 K4 = 0 0 û 0 , where det S = 1. The relations between the system variables now read T̂i(x) − P1(x) û = 0 T̂e(x) − P2(x) û = 0, (28) with P1(x) = − " (β + 4β) sinh √ 4ν2 s2+2 4ν 4β s+β 2 x √ 4ν2 s2+2 4ν 4β s+β 2 # × e− (s ν+β) x P2(x) = − " cosh √ 4ν2 s2+2 4ν 4β s+β 2 x + (s 4ν + 4β) sinh √ 4ν2 s2+2 4ν 4β s+β 2 x √ 4ν2 s2+2 4ν 4β s+β 2 # × e− (s ν+β) x . (29) We notice that the input û itself can be used as a basic output in this, parallel flow, case (Q = 1 here). In the counter-flow case we exchange β and −4β, as well as ν and −4ν, and we use vi > 0; in addition it is useful to choose the argument (x − 1). With this, we end up with the representation Q̃ T̂i(x) − P̃1(x) û = 0 Q̃ T̂e(x) − P̃2(x) û = 0 10 where Q̃ = e−(s 4ν+4β) × cosh √ ν2 s2+2 ν β s+4β2 + (s ν+β) sinh √ ν2 s2+2 ν β s+4β2 √ ν2 s2+2 ν β s+4β2 ! P̃1(x) = (β + 4β) e(s 4ν+4β) (x−1) × sinh √ ν2 s2+2 ν β s+4β2 (1−x) √ ν2 s2+2 ν β s+4β2 P̃2(x) = e(s 4ν+4β) (x−1) × cosh √ ν2 s2+2 ν β s+4β2 (1−x) + (s ν+β) sinh √ ν2 s2+2 ν β s+4β2 (1−x) √ ν2 s2+2 ν β s+4β2 ! . This system admits b̂ = Q̃−1 û as a Q̃-basic output. Using this Q̃-basic output we can express T̂i(x) and T̂e(x) via b̂ = Q̃−1 T̂e(0) without making use of the inverse of the operator Q̃: T̂i(x) = P̃1(x) b̂ T̂e(x) = P̃2(x) b̂. B. Transitions between stationary regimes One may use tables of correspondence for the Laplace transform in order to calculate the functions corresponding to the operators and operational functions occurring in the equations of the heat exchangers. For example, no. 337 in the table of G. Fodor’s book [29], tells us L−1 ( e− T √ σ2+α2 √ σ2 + α2 ) = 1 (t − T) J0 α p t2 − T2 for T > 0, where J0 is the zero order Bessel function of the first type. Extend this correspondence to T < 0 for calculating the correspondences of cosh(T √ σ2 + α2) and sinh(T √ σ2 + α2) / ( √ σ2 + α2). Substituting σ = s + 4β 4ν , α = q β 2 + (4β) 2 4ν , and T = 4ν x for the parallel flow exchanger, we end up with Ti(x, t) = (β + 4β) e− β x 2 4ν × 4ν x Z −4ν x J0 α q τ2 − (4ν x) 2 × e− 4β 4ν τ Te(0, t − νx − τ) dτ and Te(x, t) = α e− β x 2 × 4ν x Z −4ν x 4ν x−τ √ τ2−(4ν x)2 J1 (α √ τ2−(4ν x)2 ) × e− 4β 4ν τ Te(0, t − νx − τ) dτ + e(4β−β)x Te(0, t − (ν − 4ν)x), with Bessel’s function J1. This way, for both temperature fields we get integrals with bounded support, which may be called “distributed delay”. The values of Te(0, θ) for θ ∈ [t − (ν + 4ν)x), t − (ν − 4ν)x] are required to evaluate the formulae. Finally, for the counter-flow exchanger formally the same operators are involved. It is, hence, enough to replace β by −4β and ν by −4ν, and vice versa. With α = − √ β 2 +(4β)2 / ν we get Ti(x, t) = (β + 4β)e−4β(1−x) 2 ν × ν (1−x) Z −ν (1−x) J0 (α √ τ2−(ν (1−x))2 ) × e− β ν τ b(t − τ − 4ν(1 − x)) dτ Te(x, t) = αe−4β(1−x) 2 × ν (1−x) Z −ν (1−x) (ν(1−x)−τ) J1 “ α √ τ2−(ν (1−x))2 ” √ τ2−(ν(1−x))2 × e− β ν τ b(t − τ − 4ν(1 − x)) dτ + e(β−4β)(1−x) b(t + (ν − 4ν) (1 − x)). Besides a delay there is also an advance (both finite). There- fore, the trajectory of b must be known by ν in advance. For the calculation of the temperatures at x = 0 one requires b(t+ν). We observe that the flat output represents the temperature at the outflow of the exterior tube here: b(t) = Te(1, t). Remark: Analogous operators intervene in the discussion of the telegraphers equation. In [30], [31] this is exploited for the solution of a problem of signal transmission on a long electric line, where a sequence of voltage pulses at the sender’s side is computed in such a way that the influence of the line is pre-compensated. C. Trajectory planning and control In order to parameterize a transition between two stationary regimes it is now sufficient to prescribe a trajectory • [0, t∗] 3 t 7→ u(t) ∈ R to the control u for the counter- flow case; 11 • and [0, t∗] 3 t 7→ b(t) ∈ R, a trajectory for the temperature b(t) = Te(1, t) at the outflow of the exterior tube in the counter-flow case. Notice that there is no differentiability constraint on the tra- jectories (in contrast to the heat equation [4] and the examples of Sections VI and VII). Here we have a hyperbolic system. For the parallel flow exchanger the final value of the flat output is computed from the stationary interior temperature to be reached at the end of the transition, Ti(1, t∗) = Tistat. (1). For this, the stationary solution is exploited: u(t∗) = Testat. (0) = 2β (β + 4β) 1 − e− 2β Ti(1, t∗). The choice of the curve joining the two points u(0) = 0 and u(t∗) is a priori arbitrary. For the counter-flow exchanger we may proceed analogously. IX. EXAMPLE: HEAVY ROPES WITH HORIZONTALLY MOVING SUSPENSION POINT We consider the small planar deformations of heavy ropes (homogenous, with constant cross sections) with respect to a vertical rest position (in the gravity field). The suspension point can slide horizontally (such as on a gantry crane, for instance), see Fig. 5. The trajectory planning approach exposed here has been discussed by N. Petit and P. Rouchon for a single rope [32]. They have also given results for ropes carrying loads at their free end and for ropes with non-constant cross sections [32]–[34]. x x = 0 w(x, t) u(t) x = L Fig. 5. Scheme of a heavy rope. With small deformations the distance w(x, t) of the rope points w.r.t. a vertical reference satisfies the p.d.e. ∂ ∂x gx ∂w(x, t) ∂x − ∂2 w(x, t) ∂t2 = 0, x ∈ [0, L], t ≥ 0, with the boundary conditions w(L, t) = u(t), ∂w ∂x (0, t) = 0, t ≥ 0. Here, the suspension point is situated at x = L, the free end at x = 0, and u(t) is the horizontal distance of the suspension point, x being the distance of a point from the free end (we have a curved coordinate). The boundary at x = 0 corresponds to absence of a lateral force. In the beginning everything is at rest: w(x, 0) = 0, ∂w ∂t (x, 0) = 0, x ∈ [0, L]. Obviously, we have a generalization of the wave equation describing string vibrations, which results from the equation considered here by assuming constant longitudinal tension in the rope. Replacing time derivation by the operator s (zero initial conditions), we get the o.d.e. d dx gx dŵ dx (x) − s2 ŵ(x) = 0, x ∈ [0, L], and the boundary conditions ŵ(L) = û, dŵ dx (0) = 0. Transforming the independent variable with z = 2s p x/g, for x ∈ (0, L], leads to Bessel’s equation z2 d2 ŵ dz2 + z dŵ dz − z2 ŵ = 0. Thus, the general solution has the form ŵ(x) = A I0(2s p x/g) + B K0(2s p x/g), where I0 and K0 are the zero order modified Bessel functions of the first and second type, respectively. (Taking the limits enables us to extend the calculation to x = 0.) At x = 0 one has 0 = A I1(0) + B K1(0), and, as I1(0) = 0 and K1(x) tends to infinity for decreasing positive x, one also has B = 0. The boundary condition at x = L leads to û = A I0(2s p L/g), and one, thus, has I0(2s p L/g) ŵ(x) = I0(2s p x/g) û. (30) A. Introducing a Q-basic output With P(x) = I0(2s p x/g) and Q = I0(2s p L/g) = P(L), one has the system Q ŵ − P û = 0. A Q-basic output b̂ is introduced by ŵ = I0(2s p x/g) b̂ û = I0(2s p L/g) b̂. Evaluating the relation between ŵ and b̂ at x = 0 one directly gets ŵ(0) = b̂, and the Q-basic output b̂ of the system, thus, corresponds to the horizontal position of the free end of the rope. 12 B. Planning trajectories Function I0 admits a representation as a Poisson integral [35] I0(x) = 1 π Z π 2 − π 2 ex sin θ dθ, thus, I0(2s p x/g) = 1 π Z π 2 − π 2 e2s √ x/g sin θ dθ. In calculating ŵ = I0(2s p x/g)b̂ the exponential function introduces a shift (with amplitude 2 p x/g sin θ): w(x, t) = 1 π Z π 2 − π 2 b(t + 2 p x/g sin θ) dθ. As a consequence, future and past values of b are involved in the calculation, those of t − 2 p x/g and t + 2 p x/g. Here, 2 p L/g corresponds to the time required by a wave to travel from one end of the rope to the other end. As for the heat exchangers in Section VIII we have a system involving finite distributed delays and advances. Fig. 6. Positions of a heavy rope during a transition between two rest positions. C. Two ropes with a common suspension point If we consider two ropes which share a common point of suspension, generalizing (30) we immediately get the opera- tional model I0(2s p L1/g) ŵ1(x) = I0(2s p x/g) û I0(2s p L2/g) ŵ2(x) = I0(2s p x/g) û, where L1 and L2 are the lengths, w1 and w2 the distances from the initial rest position. As for the classical example of two mathematical pendulums with horizontally moving suspension point, we are confronted with a loss of controllability if the lengths are equal, L1 = L2 = L. With P(x) = I0(2s p x/g), Q = I0(2s p L/g) = P(L), we have Q ŵ1 = P û Q ŵ2 = P û and Q (ŵ1 − ŵ2) = 0 in that case. However, we can introduce b̂ in such a way that ŵ1 = ŵ2 = I0(2s p x/g) b̂ û = I0(2s p L/g) b̂. The difference of the deviations, which corresponds to the non- controllable (because autonomous) part, stays at zero during the movement. For different lengths, let be P(x) = I0(2s p x/g), x ∈ [0, L], Q1 = I0(2s p L1/g) = P(L1), and Q2 = I0(2s p L2/g) = P(L2). Introducing b̂ by ŵ1 = Q2P b̂ = I0(2s p L2/g)I0(2s p x/g) b̂ ŵ2 = Q1P b̂ = I0(2s p L1/g)I0(2s p x/g) b̂ û = Q1Q2 b̂ = I0(2s p L1/g)I0(2s p L2/g) b̂ makes it a (Q1Q2)-basic output of the system. It admits an interesting interpretation as the position of the free end of a ficticious rope of length L = L1 + L2. For the horizontal position of that ficticious rope we have ŵ(x) = I0(2s p x/g) b̂. For x = L1 ŵ(L1) = I0(2s p L1/g) b̂ = ŵ2(0), for x = L2 ŵ(L2) = I0(2s p L2/g) b̂ = ŵ1(0), and for the position of the suspension point, uL, we deduce ûL = I0(2s p L/g) b̂. Therefore, motion planning between rest positions for two ropes leads to motion planning for a single (ficticious) rope of length L = L1 + L2. If the rope positions coincide in the beginning, then the horizontal positions of the free ends of the (real) shorter ropes, ŵ2(0) and ŵ1(0), coincide with those of the points at distance L1 and L2 of the free end of the (ficticious) rope of length L1 + L2. These points are always, pairwise, on the same vertical line (and there is no other point with this property). These results, which can be easily generalized to the case of several ropes, are illustrated by a simulation result in Fig. 7. Fig. 7. Two ropes in a rest to rest motion (third ficticious rope dotted). 13 X. CONCLUSION Despite its undeniable practical importance, trajectory plan- ning for distributed parameter systems does not seem to have attracted a lot of attention in the literature before the “flatness-based”, or π-freeness, approach discussed here has been developed. After its introduction in the context of systems modeled by the linear (one dimensional) wave equation [36], [37] (see also [38], [39]), where it appeared to be most useful, it was natural to generalize the notion of π-free outputs (or π-flat outputs) to more general p.d.e. Since the first steps in this direction [9] (see also [40], [41]) quite a few cases have been studied, for example: • linear parabolic equations with constant coefficients with a single spatial variable, like the heat equation or tubular reactor models [4], [19], [20], [23], [42], • linear parabolic equations with variable coefficients with a single spatial variable or in cylindrical coordinates [15], [18], [20], [43], • Euler-Bernoulli- and Rayleigh-Bernoulli beam models [14], [25], [26], • the telegraphers equation [30], [31], • hyperbolic heat exchanger models [28], [44], • several variants of heavy rope models [32]–[34] • and the Timoshenko beam model [26], [27]. Furthermore, extension has been pursued beyond the class of linear systems discussed here, such for quasi-linear parabolic problems [20]–[22], [24], hyperbolic chemical tubu- lar reactor models [42], or the shallow water equation [45]– [47]. Other approaches generalizing flatness of finite dimen- sional nonlinear systems can be found in [34], [43], [48]–[50], for example. ACKNOWLEDGMENT This work has been partially supported by the “Deutsche Forschungsgemeinschaft”. REFERENCES [1] J. Mikusiński, Operational Calculus. Pergamon, Oxford & PWN, Warszawa, 1983, vol. 1. [2] J. Mikusiński and T. K. Boehme, Operational Calculus. Pergamon, Oxford & PWN, Warszawa, 1987, vol. 2. [3] M. Fliess, J. Lévine, P. Martin, and P. Rouchon, “Flatness and defect of non-linear systems: introductory theory and examples,” Internat. J. Control, vol. 61, pp. 1327–1361, 1995. [4] P. Martin, R. M. Murray, and P. Rouchon, “Flat systems,” in Plenary Lectures and Mini-Courses, 4th European Control Conference, Brussels, Belgium, G. Bastin and M. Gevers, Eds., 1997, pp. 211–264. [5] M. Fliess, J. Lévine, P. Martin, and P. Rouchon, “A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems,” IEEE Trans. Automat. Control, vol. AC–44, pp. 922–937, 1999. [6] H. Mounier, Propriétés structurelles des systèmes linéaires à retards : aspects théoriques et pratiques. Thèse de Doctorat, Université Paris- Sud, Orsay, 1995. [7] M. Fliess and H. Mounier, “Controllability and observability of linear delay systems: an algebraic approach,” ESAIM: Control, Optimisation and Calculus of Variations, vol. 3, pp. 301–314, 1998. [8] H. Mounier and J. Rudolph, “Flatness based control of nonlinear delay systems: A chemical reactor example,” Internat. J. Control, vol. 71, pp. 871–890, 1998. [9] M. Fliess, H. Mounier, P. Rouchon, and J. Rudolph, “Systèmes linéaires sur les opérateurs de Mikusiński et commande d’une poutre flexible,” in ESAIM Proc., vol. 2, 1997, pp. 183–193, (http://www.emath.fr/proc). [10] J. Rudolph, Beiträge zur flachheitsbasierten Folgeregelung linearer und nichtlinearer Systeme endlicher und unendlicher Dimension. Shaker Verlag, 2003. [11] ——, Flatness based control of distributed parameter systems, ser. Berichte aus der Steuerungs- und Regelungstechnik. Aachen: Shaker Verlag, 2003. [12] J. Rudolph, J. Winkler, and F. Woittennek, Flatness based control of distributed parameter systems: Examples and computer exercises from various technological domains, ser. Berichte aus der Steuerungs- und Regelungstechnik. Aachen: Shaker Verlag, 2003. [13] G. Doetsch, Handbuch der Laplace-Transformation. Basel und Stuttgart: Birkhäuser Verlag, 1971–73, vol. 1–3. [14] Y. Aoustin, M. Fliess, H. Mounier, P. Rouchon, and J. Rudolph, “Theory and practice in the motion planning and control of a flexible robot arm using Mikusiński operators,” in Proc. 5th Symposium on Robot Control, Nantes, France, 1997, pp. 287–293. [15] R. Rothfuß, U. Becker, and J. Rudolph, “Controlling a solenoid valve — a distributed parameter approach,” in Proc. 14th Int. Symp. Mathematical Theory of Networks and Systems — mtns 2000, Perpignan, France, 2000. [16] M. Gevrey, “La nature analytique des solutions des équations aux dérivées partielles,” Ann. Sci. Ecole Norm. Sup., vol. 25, pp. 129–190, 1918. [17] J.-P. Ramis, Séries divergentes et théories asymptotiques. Marseille: Soc. Math. France, 1993. [18] B. Laroche, P. Martin, and P. Rouchon, “Motion planning for the heat equation,” Int. J. Robust and Nonlinear Control, vol. 10, pp. 629–643, 2000. [19] M. Fliess, H. Mounier, P. Rouchon, and J. Rudolph, “A distributed parameter approach to the control of a tubular reactor: a multi-variable case,” in Proc. 37th IEEE Conference on Decision and Control, Tampa, FL, 1998, pp. 439–442. [20] A. F. Lynch and J. Rudolph, “Flachheitsbasierte Randsteuerung parabolischer Systeme mit verteilten Parametern,” at — Automa- tisierungstechnik, vol. 48, pp. 478–486, 2000. [21] ——, “Flatness-based boundary control of a nonlinear parabolic equa- tion modelling a tubular reactor,” in Nonlinear Control in the Year 2000 (Vol. 2), ser. Lecture Notes in Control and Inform. Sci., A. Isidori, F. Lamnabhi-Lagarrique, and W. Respondek, Eds. Springer-Verlag, 2000, vol. 259, pp. 45–54. [22] ——, “Flatness-based boundary control of a class of quasilinear parabolic distributed parameter systems,” Internat. J. Control, 2002. [23] M. Fliess, H. Mounier, P. Rouchon, and J. Rudolph, “Controlling the transient of a chemical reactor: a distributed parameter approach,” in Proc. Computational Engineering in Systems Application IMACS Multiconference, (CESA’98), Hammamet, Tunisia, 1998. [24] A. F. Lynch and J. Rudolph, “Flatness-based boundary control of two coupled nonlinear pdes modelling a tubular reactor,” in Proc. of the International Symposium on Nonlinear Theory and its Applications (NOLTA 2000), vol. 2, Dresden, Germany, 2000, pp. 641–644. [25] W. Haas and J. Rudolph, “Steering the deflection of a piezoelectric ben- der,” in Proc. 5th European Control Conference, Karlsruhe, Germany, 1999. [26] J. Rudolph and F. Woittennek, “Flachheitsbasierte Randsteuerung von elastischen Balken mit Piezoaktuatoren,” at — Automatisierungstechnik, vol. 50, pp. 412–421, 2002. [27] ——, “Flachheitsbasierte Steuerung eines Timoshenko-Balkens,” Z. angew. Math. Mech., vol. 83, pp. 119–127, 2003. [28] J. Rudolph, “Commande par platitude de modèles à paramètres répartis pour échangeurs de chaleur,” in Actes CIFA 2000, Ecole Centrale de Lille, 5–8 juillet 2000, 2000, pp. 566–571. [29] G. Fodor, Laplace Transforms in Engineering. Budapest: Akadémiai Kiadó, 1965. [30] M. Fliess, P. Martin, N. Petit, and P. Rouchon, “Commande de l’équation des télégraphistes et restauration active d’un signal,” Traitement du Signal, vol. 15, Spécial, pp. 619–625, 1998. [31] ——, “Active signal restoration for the telegraph equation,” in Proc. 38th IEEE Conference on Decision and Control, 1999. [32] N. Petit and P. Rouchon, “Motion planning for heavy chain systems,” in Nonlinear Control in the Year 2000 (Vol. 2), ser. Lecture Notes in Control and Inform. Sci., A. Isidori, F. Lamnabhi-Lagarrique, and W. Respondek, Eds. Springer-Verlag, 2000, vol. 259, pp. 229–236. [33] ——, “Motion planning for heavy chain systems,” SIAM J. Control Optim., vol. 40, pp. 475–495, 2001. [34] N. Petit, Systèmes à retards. Platitude en génie des procédés et contrôle de certaines équations des ondes. Thèse de Doctorat, École Nationale Supérieure des Mines de Paris, 2000. 14 [35] G. N. Watson, A Treatise on the Theory of Bessel Functions. Cambridge University Press, 1922 (Reprint 1996). [36] H. Mounier, J. Rudolph, M. Petitot, and M. Fliess, “A flexible rod as a linear delay system,” in Proc. 3rd European Control Conference, Rome, Italy, 1995, pp. 3676–3681. [37] M. Fliess, H. Mounier, P. Rouchon, and J. Rudolph, “Controllability and motion planning for linear delay systems with an application to a flexible rod,” in Proc. 34th IEEE Conference on Decision and Control, New Orleans, 1995, pp. 2046–2051. [38] H. Mounier, P. Rouchon, and J. Rudolph, “Some examples of linear systems with delays,” RAIRO—APII—JESA, vol. 31, pp. 911–925, 1997. [39] H. Mounier, J. Rudolph, M. Fliess, and P. Rouchon, “Tracking control of a vibrating string with an interior mass viewed as delay system,” ESAIM: Control, Optimisation and Calculus of Variations, vol. 3, pp. 315–321, 1998. [40] M. Fliess and H. Mounier, “Tracking control and π-freeness of infinite dimensional linear systems,” in Dynamical Systems, Control, Coding, Computer Vision, G. Picci and D. S. Gilliam, Eds. Birkhäuser, 1999, pp. 45–68. [41] M. Fliess, “Variations sur la notion de contrôlabilité,” in Quelques Aspects de la Théorie du Contrôle. Journée Annuelle Soc. Math. France, Paris, France, 2000, pp. 47–86. [42] P. Rouchon and J. Rudolph, “Réacteurs chimiques différentiellement plats : planification et suivi de trajectoires,” in Commande de procédés chimiques — Réacteurs et colonnes de distillation, J. P. Corriou, Ed. Hermès Science Publications, 2001, ch. 5, pp. 163–200. [43] B. Laroche, Extension de la notion de platitude à des systèmes décrits par des équations aux dérivées partielles linéaires. Thèse de Doctorat, École Nationale Supérieure des Mines de Paris, 2000. [44] J. Rudolph, “Randsteuerung von Wärmetauschern mit örtlich verteilten Parametern: Ein flachheitsbasierter Zugang,” at — Automatisierungs- technik, vol. 48, pp. 399–406, 2000. [45] F. Dubois, N. Petit, and P. Rouchon, “Motion planning and nonlinear simulations for a tank containing a fluid,” in Proc. 5th European Control Conference, Karlsruhe, Germany, 1999. [46] P. Rouchon, “Dynamique et contrôle des systèmes : stabilité, équivalence et systèmes plats,” Mémoire présenté à l’Université de Paris-Sud pour obtenir le diplôme d’habilitation à diriger des recherches, septembre 2000. [47] N. Petit and P. Rouchon, “Dynamics and solutions to some control problems for water-tank systems,” IEEE Trans. Automat. Control, vol. AC–47, pp. 594–609, 2002. [48] P. Martin, R. M. Murray, and P. Rouchon, “Flat systems: open problems, infinite dimensional extension, symmetries and catalog,” in Advances in the Control of Nonlinear Systems, ser. Lecture Notes in Control and Inform. Sci., A. Baños, F. Lamnabhi-Lagarrigue, and F. J. Montoya, Eds. Springer-Verlag, 2001, vol. 264, pp. 33–57. [49] F. Ollivier and A. Sedoglavic, “A generalization of flatness to nonlinear systems of partial differential equations. Applications to the command of a flexible rod,” in Proc. Fifth IFAC Symposium on Nonlinear Control Systems (NOLCOS 2001), Saint-Petersburg, Russia, 2001. [50] P. Rouchon, “Motion planning, equivalence, infinite dimensional sys- tems,” Int. J. Appl. Math. Comput. Sci., vol. 11, pp. 165–188, 2001. PLACE PHOTO HERE Joachim Rudolph Dr. Rudolph got an engineering diploma in technical cybernetics from the “Uni- versität Stuttgart”, Germany in 1989 and a “Doc- torat” from the “Université Paris XI, Orsay”, France in 1991. He obtained the “Dr.-Ing. habil.” degree and the title “Privatdozent” from “Technische Uni- versität Dresden”, Germany in 2003. He is a se- nior researcher (“Oberassistent”) at the “Institut für Regelungs- und Steuerungstheorie” of that univer- sity. His main research interests are in controller and observer design for nonlinear systems, algebraic sys- tems theory, and infinite dimensional systems. Practical applications presently concern the use of magnetic bearings in high precision tooling spindles and crystal growth processes.