Backstepping Control of Coupled Linear Parabolic PIDEs with Spatially-Varying Coefficients
IIEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 1
Backstepping Control of Coupled Linear ParabolicPIDEs with Spatially-Varying Coefficients
Joachim Deutscher and Simon Kerschbaum
Abstract —This paper considers the backstepping design ofstate feedback controllers for coupled linear parabolic par-tial integro-differential equations (PIDEs) of Volterra-type withdistinct diffusion coefficients, spatially-varying parameters andmixed boundary conditions. The corresponding target systemis a cascade of parabolic PDEs with local couplings allowinga direct specification of the closed-loop stability margin. Thedetermination of the state feedback controller leads to kernelequations, which are a system of coupled linear second-orderhyperbolic PIDEs with spatially-varying coefficients and ratherunusual boundary conditions. By extending the method of succes-sive approximations for the scalar case to the considered systemclass, the well-posedness of these kernel equations is verifiedby providing a constructive solution procedure. This results ina systematic method for the backstepping control of coupledparabolic PIDEs as well as PDEs. The applicability of the newbackstepping design method is confirmed by the stabilization oftwo coupled parabolic PIDEs with Dirichlet/Robin unactuatedboundaries and a coupled Neumann/Dirichlet actuation.
Index Terms —Distributed-parameter systems, parabolic sys-tems, coupled PIDEs, boundary control, backstepping.
I. I
NTRODUCTION I N the last decade the backstepping approach has emergedas a very powerful tool to solve stabilization problems forboundary controlled distributed-parameter systems (DPS) (see[12] and [15] for an overview). The basic idea of this method isto utilize an invertible Volterra-type integral transformation forfacilitating the controller design. As far as parabolic systemsare concerned, the backstepping method for the state feedbackdesign was first introduced in [14], [21] for scalar parabolicPDEs and PIDEs. Subsequently, these results were extended toparabolic systems with space- and time-dependent coefficientsin [22], [16], to parabolic PDEs with Volterra nonlinearities in[24], [25] and to higher-dimensional spatial domains in [15],[9].Recently, the backstepping control of coupled parabolicsystems attracted the attention of many researchers. This systemclass is not only of paramount theoretical appeal, but it alsohas a major importance for applications. Typical real-worldproblems originate from chemical and biochemical engineering,whereby chemical fixed-bed and tubular reactors are importantexamples (see, e. g., [20], [10]). A first solution to the backstep-ping design for coupled diffusion-reaction systems with equaldiffusion coefficients and constant coefficients was given in [1],[2], [3]. In this work it was shown how the approach of [21] fordetermining the scalar integral kernel, that defines the backstep-ping transformation, can directly be extended to the matrix case.
J. Deutscher and S. Kerschbaum are with the Lehrstuhl f¨ur Regelungstechnik,Universit¨at Erlangen-N¨urnberg, Cauerstraße 7, D-91058 Erlangen, Germany.(e-mail: { joachim.deutscher,simon.kerschbaum } @fau.de) The latter appears in the control of n coupled PDEs, becausethe state vector is n -dimensional. More challenging is thebackstepping control of coupled parabolic systems with distinctdiffusion coefficients. First solutions to this problem considereddiffusion-reaction systems with constant coefficients (see [3],[13], [19]). For this system class it is possible to simplifythe kernel equations by assuming a diagonal structure of thematrix kernel. However, this approach cannot be extended tothe spatially-varying coefficients case. Recently, a very generaland elegant solution of the backstepping problem for diffusion-convection-reaction systems with spatially-varying coefficientsand Dirichlet boundary conditions (BCs) was presented in[26]. Thereby, no assumption on the form of integral kernel isimposed. By introducing a local coupling term with a triangularstructure in the target system it was shown that the resultingkernel equations are well-posed. The latter result followed fromestablishing an interesting relation of the kernel equations forthe considered system class with the kernel equations relatedto coupled first-order hyperbolic PDEs treated in [7], [8].A larger class of DPS can be modelled by partial integro-differential equations (PIDEs) . They arise directly in thephysical modelling or from a singular perturbation of separatesubsystems with different time-scales. For this system class,Robin BCs have to be taken into account frequently. Thisresults, e. g., from the energy balance along the boundariesin heat transfer problems if the DPS is in direct contactwith a surrounding medium. Accordingly, heat isolation at theboundaries requires to introduce Neumann BCs. Moreover, thelinearisation of a nonlinear DPS around a steady state, non-homogenous materials or unusually shaped spatial domainslead to spatially-varying coefficients. Hence, it is of interestto extend the backstepping method to this class of coupledparabolic systems.This paper is concerned with the backstepping control ofparabolic PIDEs of Volterra-type with distinct diffusion coef-ficients and spatially-varying parameters. Thereby, Dirichlet,Neumann or Robin BCs are assumed independently for eachcomponent of the state. Similar to [26] a target system witha local coupling term in triangular form is proposed. Thisresults in a cascaded set of parabolic PDEs, which allow asimple proof of well-posedness and stability. Moreover, thecorresponding rate of exponential convergence can directly bespecified in the design so that a simple parametrization ofthe target system is possible. The main contribution of thepaper is the verification of the well-posedness of the resulting kernel equations . For n coupled parabolic PIDEs they are asystem of n coupled hyperbolic PIDEs with a Klein-Gordon-type spatial differential operator, spatially-varying coefficientsand rather unusual boundary conditions with integral terms. a r X i v : . [ m a t h . O C ] D ec EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 2
Since this type of kernel equations cannot be traced backto the kernel equations found for general heterodirectionalhyperbolic systems in [7], [8] with the result in [26], a newmethod for their solution is needed. In order to provide asystematic solution procedure, the approach of [21], [22] isdirectly extended to the system class in question. The resultingkernel equations for the diagonal elements share the same formas in the scalar case. However, the off-diagonal elements aregoverned by kernel equations with a different structure. In afirst step, the kernel PIDEs with spatially-varying coefficientsare simplified utilizing suitable transformations. In particular,by extending the results in [22], new kernel PIDEs are obtained,where the coefficients w. r. t. the second-order partial derivativesare equal to one. Hence, the usual linear change of coordinatesmapping the kernel PIDE into its canonical form can directly beapplied. The resulting kernel equations allow a very systematicformulation of the corresponding integral equations. This leadsto the usual result with the standard triangular spatial domainfor the diagonal kernel elements. However, the spatial domainrelated to the off-diagonal kernel PIDEs is no longer restrictedto the first quadrant. It is shown that by introducing suitableartificial BCs, well-posed kernel equations are obtained. Thelatter property is proved by verifying the uniform convergenceof the corresponding successive approximation. This resultsin a systematic method for the backstepping state feedbackstabilization of the coupled PIDEs in question.As far as coupled systems of parabolic PDEs are concerned,the presented approach provides a direct solution of the n second-order hyperbolic kernel PDEs for Dirichlet, Neumannand Robin BCs. In contrast, the recent approach in [26]assuming Dirichlet BCs solves a system of n first-ordertransport equations, which follow from the second-order kernelPDEs. The new approach could also be of interest whenconsidering other types of coupled PDEs. For example, itmay be the immediate starting point for an extension tocoupled parabolic systems with spatially- and temporally-varying reaction by directly utilizing the results [22], [16] forscalar parabolic systems.The next section formulates the considered problem. Then,the target system is introduced and the corresponding well-posedness and stability is proved in the following section. Asystematic procedure for the solution of the kernel equations ispresented in Section IV. Subsequently, Section V investigatesthe stability of the resulting closed-loop system. The proposedbackstepping-based state feedback controller design is illus-trated for an unstable coupled parabolic system of two PIDEs,where a Dirichlet and Robin BC is imposed at the unactuatedboundary and the actuation is of Neumann and Dirichlet typewith a coupling at that boundary. II. P ROBLEM FORMULATION
Consider the following system described by coupled linearparabolic PIDEs ∂ t x ( z, t ) = Λ( z ) ∂ z x ( z, t ) + Φ( z ) ∂ z x ( z, t ) + A ( z ) x ( z, t )+ A ( z ) x (0 , t ) + (cid:90) z F ( z, ζ ) x ( ζ, t )d ζ (1a) θ [ x ( t )](0) = 0 , t > (1b) θ [ x ( t )](1) = u ( t ) , t > (1c)with (1a) defined on ( z, t ) ∈ (0 , × R + , the state x ( z, t ) ∈ R n , n > , and the input u ( t ) ∈ R n . The matrix Λ( z ) ∈ R n × n given by Λ( z ) = diag( λ ( z ) , . . . , λ n ( z )) (2)contains mutually different, positive and spatially-varying diffusion coefficients λ i ∈ C [0 , , i = 1 , , . . . , n , i. e., λ i ( z ) (cid:54) = λ j ( z ) , z ∈ [0 , , i (cid:54) = j . The convection term ischaracterized by the diagonal matrix Φ( z ) = diag(Φ ( z ) , . . . , Φ n ( z )) (3)with Φ ∈ ( C [0 , n × n and the matrix A = [ A ij ] ∈ ( C [0 , n × n describes the reaction term . The local term inthe PIDE (1a) is determined by A ∈ ( C [0 , n × n and the non-local term by the integral kernel F ∈ ( C ([0 , )) n × n .The BCs are represented by the formal matrix differentialoperators θ i [ h ] = B i d z h + B i h, i = 0 , . (4)Therein, B = (cid:20) m I p (cid:21) and B = (cid:20) I m Q (cid:21) , (5) m + p = n , is assumed, which means that the first m BCs at z = 0 are of Dirichlet type. Accordingly, the remaining p BCsrepresent Neumann/Robin BCs, in which Q = diag( q , . . . , q p ) ∈ R p × p (6)is a diagonal matrix. This setup can always be achieved fordecoupled BCs at z = 0 by a suitable reordering of the state x . Moreover, the actuation at z = 1 is described by a diagonalmatrix B ∈ R n × n and an arbitrary matrix B ∈ R n × n . Thisspecifies Dirichlet, Neumann or Robin actuation at z = 1 independently for each component u i , i = 1 , , . . . , n , ofthe input u . Thereby, also different types of BCs for eachcomponent x i ( z, t ) ∈ R , i = 1 , , . . . , n , of the state x at z = 1 as well as possible couplings of these BCs are included. Remark 1:
The type of convection and BCs specified by(3) and (4) appear frequently in applications. Well-knownexamples are chemical fixed-bed and tubular reactors (see,e. g., [20], [10]). (cid:47)
The initial conditions (ICs) of the system are x ( z,
0) = x ( z ) with x ∈ ( L (0 , n .The system (1) can be simplified by introducing the bound-edly invertible Hopf-Cole-type state transformation ˇ x ( z, t ) = exp (cid:16) (cid:90) z Λ − ( ζ )Φ( ζ )d ζ (cid:17) x ( z, t ) . (7) EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 3
With a direct calculation it is easy to verify that (7) maps (1a)into the coupled PIDEs ∂ t ˇ x ( z, t ) = Λ( z ) ∂ z ˇ x ( z, t ) + ˇ A ( z )ˇ x ( z, t )+ ˇ A ( z )ˇ x (0 , t ) + (cid:90) z ˇ F ( z, ζ )ˇ x ( ζ, t )d ζ (8)on ( z, t ) ∈ (0 , × R + . Thereby, the resulting system sharesthe same type of BCs as (1). Consequently, the convection termin (1) is omitted in the sequel, i. e., Φ( z ) ≡ is assumed.Consider the state feedback controller u ( t ) = K [ x ( t )] (9)with the formal feedback operator K mapping the state x ( z, t ) ∈ R n to the input u ( t ) ∈ R n . The problem consideredin this paper is the backstepping design of (9) such that theresulting closed-loop system is exponentially stable with aprescribed rate of convergence.III. S TATE FEEDBACK CONTROLLER DESIGN
A. Selection of the Target System
In what follows, the backstepping approach (see, e. g., [12])is used to determine the state feedback controller (9). To thisend, the backstepping coordinates ˜ x ( z, t ) = x ( z, t ) − (cid:90) z K ( z, ζ ) x ( ζ, t )d ζ (10)with the integral kernel K ( z, ζ ) ∈ R n × n are introduced for (1).The determination of the feedback gains in (9) is significantlysimplified in these coordinates. As a first step in the design,an exponentially stable target system has to be found so thatthe change of coordinates (10) exists. By following the linesof [21] for the scalar case and the recent results in [7], [26],the following target system ∂ t ˜ x ( z, t ) = Λ( z ) ∂ z ˜ x ( z, t ) − µ c ˜ x ( z, t ) (11a) − (cid:101) A ( z )( E E (cid:62) ∂ z ˜ x (0 , t ) + E E (cid:62) ˜ x (0 , t )) θ [˜ x ( t )](0) = 0 , t > (11b) ˜ θ [˜ x ( t )](1) = 0 , t > (11c)with µ c ∈ R is proposed. Therein, (11a) is defined on ( z, t ) ∈ (0 , × R + and (cid:101) A ( z ) = [ (cid:101) A ,ij ( z )] ∈ R n × n is amatrix containing only n ( n − / non-zero elements, i. e., (cid:101) A ,ij ( z ) = (cid:40) f ij ( z ) , λ i < λ j , else (12)with f ij ( z ) determined by the kernel K ( z, ζ ) . Here, λ i < λ j is the shorthand notation for λ i ( z ) < λ j ( z ) , z ∈ [0 , .The non-zero elements (cid:101) A ,ij ( z ) , λ i < λ j , are determinedby the solution of the kernel equations (cf. Section IV-A).Furthermore, the matrices E = (cid:20) I m (cid:21) ∈ R n × m and E = (cid:20) I p (cid:21) ∈ R n × p (13)were used in (11a). The BCs of the target system at z = 1 arecharacterized by ˜ θ [ h ] = (cid:101) B d z h + (cid:101) B h with (cid:101) B = diag(˜ b , . . . , ˜ b n ) and (cid:101) B = diag(˜ b , . . . , ˜ b n ) (14) being diagonal matrices satisfying | ˜ b i | + | ˜ b i | > , i =1 , , . . . , n . This means that decoupled BCs at z = 1 areassigned in the target system. Thereby, Dirichlet BCs areassigned for z = 1 , if the plant has this type of BC for theconsidered element of x . On the contrary, Neumann or RobinBCs can be interchanged. As there are p Neumann/Robin BCsin (1b), the local coupling term in (11a) has to depend on ˜ x (0 , t ) . In contrast, for Dirichlet BCs it was shown in [26]that the corresponding term requires a dependency on ∂ z ˜ x (0 , t ) .These couplings have to be introduced, because a completedecoupling of the PDEs in (11a) leads to an overdeterminedand thus unsolvable set of kernel equations.For the stability analysis of the target system (11) it isconvenient to introduce the new state ˜ x ∗ ( z, t ) = P ˜ x ( z, t ) . (15)Therein, P ∈ R n × n is a permutation matrix with P − = P (cid:62) such that the relation λ ∗ ( z ) > . . . > λ ∗ n ( z ) holds forthe diagonal elements in Λ ∗ ( z ) = diag( λ ∗ ( z ) , . . . , λ ∗ n ( z )) = P Λ( z ) P (cid:62) . Then, the target system (11) takes the form ∂ t ˜ x ∗ ( z, t ) = Λ ∗ ( z ) ∂ z ˜ x ∗ ( z, t ) − µ c ˜ x ∗ ( z, t ) (16a) − (cid:101) A ∗ ( z )( R ∂ z ˜ x ∗ (0 , t ) + R ˜ x ∗ (0 , t )) θ [ P (cid:62) ˜ x ∗ ( t )](0) = 0 (16b) ˜ θ [ P (cid:62) ˜ x ∗ ( t )](1) = 0 (16c)with R = P E E (cid:62) P (cid:62) and R = P E E (cid:62) P (cid:62) . Due tothe previous definition (12) of (cid:101) A ( z ) and the introducedpermutation of the state, the matrix (cid:101) A ∗ ( z ) is strictly lowertriangular, i. e., (cid:101) A ∗ ( z ) = P (cid:101) A ( z ) P (cid:62) = . . . . . . (cid:101) A ∗ , ( z ) . . . . . . ...... . . . . . . ... (cid:101) A ∗ ,n ( z ) . . . (cid:101) A ∗ ,n n − ( z ) 0 . (17)This implies that (16) has a cascade structure, which signifi-cantly simplifies the corresponding stability analysis. The nexttheorem shows that the target system (16) and thus (11) isexponentially stable with a stability margin assignable by µ c . Theorem 1 (Stability of the target system):
Assume that µ c > µ max , in which µ max is the largest eigenvalueof (16) for µ c = 0 and (cid:101) A ∗ ( z ) ≡ . Then, the targetsystem (16) is exponentially stable in the L -norm (cid:107) h (cid:107) = ( (cid:82) (cid:107) h ( z ) (cid:107) C n d z ) / , i. e., (cid:107) ˜ x ( t ) (cid:107) ≤ (cid:102) M e ( µ max − µ c ) t (cid:107) ˜ x (0) (cid:107) , t ≥ (18)for all ˜ x (0) ∈ ( H (0 , n compatible with the BCs (16b),(16c) and an (cid:102) M ≥ .The proof of this result can be found in the appendix. B. Derivation of the Kernel Equations
The equations to be solved for determining K ( z, ζ ) in (10)result from requiring that (10) and a suitable feedback (9) map(1) into the target system (11). Differentiating (10) w. r. t. time, EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 4 inserting (1a), utilizing (11a) and interchanging the order ofintegration in the double integral results in ∂ t ˜ x ( z, t ) = Λ( z ) ∂ z ˜ x ( z, t ) − µ c ˜ x ( z, t ) − (cid:101) A ( z )( E E (cid:62) ∂ z ˜ x (0 , t ) + E E (cid:62) ˜ x (0 , t ))+ Λ( z ) ∂ z (cid:90) z K ( z, ζ ) x ( ζ, t )d ζ + (cid:16) (cid:101) A ( z ) E E (cid:62) + A ( z ) − (cid:90) z K ( z, ζ ) A ( ζ )d ζ (cid:17) x (0 , t )+ (cid:101) A ( z ) E E (cid:62) ∂ z x (0 , t ) + ( A ( z ) + µ c I ) x ( z, t ) − (cid:90) z K ( z, ζ )Λ( ζ ) ∂ ζ x ( ζ, t )d ζ − (cid:90) z (cid:16) K ( z, ζ )( A ( ζ ) + µ c I )+ (cid:90) zζ K ( z, ¯ ζ ) F (¯ ζ, ζ )d¯ ζ (cid:17) x ( ζ, t )d ζ, (19)in which K (0 ,
0) = 0 was assumed. It is convenient to definethe integral operators B [ K ]( z, ζ ) = K ( z, ζ )( A ( ζ ) + µ c I ) − F ( z, ζ ) + (cid:90) zζ K ( z, ¯ ζ ) F (¯ ζ, ζ )d¯ ζ (20)and C [ K ]( z ) = A ( z ) − (cid:90) z K ( z, ζ ) A ( ζ )d ζ. (21)Then, integration by parts and using the Leibniz differentiationrule yield after simple calculations ∂ t ˜ x ( z, t ) = Λ( z ) ∂ z ˜ x ( z, t ) − µ c ˜ x ( z, t ) − (cid:101) A ( z ) (cid:0) E E (cid:62) ∂ z ˜ x (0 , t ) + E E (cid:62) ˜ x (0 , t ) (cid:1) + (cid:90) z (cid:16) Λ( z ) ∂ z K ( z, ζ ) − ∂ ζ ( K ( z, ζ )Λ( ζ )) − B [ K ]( z, ζ ) (cid:17) x ( ζ, t )d ζ + M ( z ) x (0 , t ) + M ( z ) ∂ z x (0 , t )+ (cid:16) Λ( z ) K (cid:48) ( z, z ) + Λ( z ) ∂ z K ( z, z )+ ∂ ζ K ( z, z )Λ( z ) + K ( z, z )Λ (cid:48) ( z )+ A ( z ) + µ c I (cid:17) x ( z, t )+ (cid:0) Λ( z ) K ( z, z ) − K ( z, z )Λ( z ) (cid:1) ∂ z x ( z, t ) (22)with ( · ) (cid:48) = dd z ( · ) . Therein, M ( z ) = (cid:101) A ( z ) E E (cid:62) + C [ K ]( z ) − ∂ ζ K ( z, − K ( z, (cid:48) (0) and M ( z ) = (cid:101) A ( z ) E E (cid:62) + K ( z, were utilized. In order to derivethe kernel BCs use E E (cid:62) + E E (cid:62) = I n and x i = E (cid:62) i x , i = 1 , , as well as the condition M ( z ) x (0 , t ) + M ( z ) ∂ z x (0 , t )= M ( z ) E x (0 , t ) + M ( z ) E x (0 , t )+ M ( z ) E ∂ z x (0 , t ) + M ( z ) E ∂ z x (0 , t ) = 0 (23)implied by the requirement that (22) coincides with (11a).Observe that (1b) gives x (0 , t ) = 0 and ∂ z x (0 , t ) = − Q x (0 , t ) in view of (4) and (5). After inserting this in (23) and comparing the BC (1b) with (11b), the kernel equations Λ( z ) ∂ z K ( z, ζ ) − ∂ ζ ( K ( z, ζ )Λ( ζ )) = B [ K ]( z, ζ ) (24a) K ( z, E = − (cid:101) A ( z ) E (24b)(24c) = (cid:0) (cid:101) A ( z )+ C [ K ]( z ) (cid:1) E + K ( z, z )Λ (cid:48) ( z ) = − ( A ( z )+ µ c I ) (24d) K ( z, z )Λ( z ) − Λ( z ) K ( z, z ) = 0 (24e) K (0 ,
0) = 0 ∂ ζ K ( z, E + K ( z, (cid:48) (0) E + Λ(0) E Q )Λ( z ) K (cid:48) ( z, z )+Λ( z ) ∂ z K ( z, z )+ ∂ ζ K ( z, z )Λ( z ) (24f)are obtained, in which (24a) is defined on < ζ < z < . Remark 2:
If all BCs at z = 0 are of the same type, thenthe substitution E → I n and E → for Dirichlet BCsand E → and E → I n for Neumann/Robin BCs in (24)yields the related kernel equations. With this, their solutionalso follows from the subsequent results. (cid:47) If the kernel K ( z, ζ ) is known, then the feedback controller(9) follows from imposing the BC (11c). For this, consider(1c) with (4) in the form u ( t ) = θ [ x ( t )](1) = B ∂ z x (1 , t ) + B x (1 , t ) . (25)If an element in the coefficient matrix (cid:101) B of ˜ θ [˜ x ( t )](1) (see(11c)) is zero, i. e., ˜ b i = 0 , then ˜ x i (1 , t ) = 0 and (10) yieldsfor the i -th component of the state vector x i (1 , t ) = (cid:90) e (cid:62) i K (1 , z ) x ( z, t )d z. (26)Therein, e i ∈ R n denotes the i -th unit vector. In the case ˜ b i (cid:54) = 0 , one has ∂ z ˜ x i (1 , t ) = − (˜ b i / ˜ b i )˜ x i (1 , t ) and obtains ∂ z x i (1 , t ) = ( e (cid:62) i K (1 , − ˜ b i ˜ b i e (cid:62) i ) x (1 , t )+ (cid:90) e (cid:62) i ( ∂ z K (1 , z ) + ˜ b i ˜ b i K (1 , z )) x ( z, t )d z (27)after differentiating (10) w. r. t. z . With this, the feedbackoperator K in (9) follows from inserting (26) and (27) in (25).Thereby, all spatial derivatives of the state possibly appearingin (25) are removed by the substitution (27).In the next section the solvability of (24) is investigated. Tothis end, (24) is converted into integral equations which canbe solved by a fixpoint iteration. By considering a sufficientlylarge but finite number of iterations this leads to the methodof successive approximations , which results in a systematicapproach for determining the kernel K ( z, ζ ) .IV. S OLUTION OF THE KERNEL EQUATIONS
A. Component Form of the Kernel Equations
For converting (24) into integral equations the boundaryvalue problems (BVPs) for the elements K ij ( z, ζ ) of the kernel K ( z, ζ ) = [ K ij ( z, ζ )] , i, j = 1 , , . . . , n , satisfying (24) arederived. Thereby, B = [ B ij ] and C = [ C ij ] are utilized. Thisyields the kernel PIDEs λ i ( z ) ∂ z K ij ( z, ζ ) − ∂ ζ ( λ j ( ζ ) K ij ( z, ζ )) = B ij [ K ]( z, ζ ) . (28) EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 5
Evaluating (24b) for λ i ≥ λ j and j ≤ m leads to the BC K ij ( z,
0) = 0 (29)in view of (12). Similarly, for λ i ≥ λ j and j > m , (24c) with(6) results in the BC λ j (0) ∂ ζ K ij ( z,
0) + ( λ (cid:48) j (0) + q j λ j (0)) K ij ( z,
0) = C ij [ K ]( z ) , (30)whereas the relations (cid:101) A ,ij ( z ) = − λ j (0) K ij ( z, , j ≤ mλ j (0) ∂ ζ K ij ( z, − C ij [ K ]( z )+ (cid:0) λ (cid:48) j (0) + q j λ j (0) (cid:1) K ij ( z, , j > m (31a)(31b)obtained for λ i < λ j determine the elements of (cid:101) A ( z ) in (12).Next, considering (24e) for i = j results in ( λ i ( z ) − λ i ( z )) K ii ( z, z ) = 0 , (32)which shows that no BC has to be fulfilled for K ii ( z, z ) tosatisfy (24e). In contrast, the case i (cid:54) = j gives ( λ j ( z ) − λ i ( z )) K ij ( z, z ) = 0 (33)so that for unequal diffusion coefficients λ i ( z ) and λ j ( z ) theBC K ij ( z, z ) = 0 (34)follows. After differentiating (34) one obtains K (cid:48) ij ( z, z ) = 0 (35)and thus ∂ ζ K ij ( z, z ) = − ∂ z K ij ( z, z ) . (36)Finally, the remaining BC (24d) for i = j leads to the ODE λ i ( z ) K (cid:48) ii ( z, z ) + λ (cid:48) i ( z ) K ii ( z, z ) = − ( A ii ( z ) + µ c ) . (37)With the IC (24f), the corresponding solution is K ii ( z, z ) = − (cid:112) λ i ( z ) (cid:90) z A ii ( ζ ) + µ c (cid:112) λ i ( ζ ) d ζ. (38)If i (cid:54) = j , then (24d) gives ( λ i ( z ) − λ j ( z )) ∂ z K ij ( z, z ) = − A ij ( z ) (39)in view of (35) and (36). Hence, the BC ∂ z K ij ( z, z ) = A ij ( z ) λ j ( z ) − λ i ( z ) (40)follows. By collecting the preceding results, the followingkernel equations for the elements of K ( z, ζ ) can be deduced.Thereby, the BC in [ · ] ∗ has to be considered only for thoseelements of the corresponding expression, for which thecondition ∗ is fulfilled. With this, the component form i = j : λ i ( z ) ∂ z K ii ( z, ζ ) − ∂ ζ ( λ i ( ζ ) K ii ( z, ζ )) = B ii [ K ]( z, ζ ) (41a) [ K ii ( z,
0) = 0] i ≤ m (41b) [ λ i (0) ∂ ζ K ii ( z, λ (cid:48) i (0)+ q i λ i (0)) K ii ( z, C ii [ K ]( z )] i>m (41c) K ii ( z, z ) = − (cid:112) λ i ( z ) (cid:90) z A ii ( ζ ) + µ c (cid:112) λ i ( ζ ) d ζ (41d) and i (cid:54) = j : λ i ( z ) ∂ z K ij ( z, ζ ) − ∂ ζ ( λ j ( ζ ) K ij ( z, ζ )) = B ij [ K ]( z, ζ ) (42a) [ K ij ( z,
0) = 0] λ i >λ j , j ≤ m (42b) [ λ j (0) ∂ ζ K ij ( z,
0) + ( λ (cid:48) j (0) + q j λ j (0)) K ij ( z, C ij [ K ]( z )] λ i >λ j , j>m (42c) K ij ( z, z ) = 0 (42d) ∂ z K ij ( z, z ) = A ij ( z ) λ j ( z ) − λ i ( z ) (42e)of the kernel equations (24) is obtained, whereby (41a) and(42a) are defined on < ζ < z < . Remark 3:
It is interesting to note that the kernel equationsfor the diagonal elements (41) consist of the usual but coupledkernel equations found in the scalar case (see [21], [22]). Incontrast, the BVP (42) for all other elements are different w. r. t.the kernel PIDE and the corresponding BCs, which needs anew approach for their solution. This is the topic of the nextsections. (cid:47)
Remark 4:
In the case of equal diffusion coefficients , i. e., λ max ≥ λ ( z ) = λ ( z ) = . . . = λ n ( z ) ≥ λ min > , z ∈ [0 , ,the kernel equations (42) have the same form as (41) and thusare easier to solve. In particular, (33) leads to no conditionfor K ij ( z, z ) , i. e., (42d) does not appear. Then, (24d) yieldsan ODE for K ij ( z, z ) replacing (42e) by a BC of the form(41d). As a consequence, (29) and (30) can be imposed for all i (cid:54) = j as kernel BCs so that (cid:101) A ( z ) ≡ holds in (11a). Thisresult was first shown in [1] for a diffusion-reaction systemwith constant coefficients and Neumann BCs. (cid:47) B. Transformation of the Kernel Equations
In what follows, the transformation approach in [22] for asingle kernel PIDE is extended to the coupled system of kernelPIDEs (41a) and (42a). Then, the well-posedness of the BVPsrelated to the resulting simpler kernel PIDEs is much easierto prove.For the ease of presentation, matrix elements and coordinateswith indices i and j may be represented without the index butbold face, e. g., KKK = K ij , holds in the sequel. Moreover, indexnotation for derivatives is used, i. e., f z = ∂ z f .
1) Transformation of the second-order partial derivatives:
In order to eliminate the dependency on z and ζ of thecoefficients w. r. t. the second-order partial derivatives in (41a)and (42a), introduce the change of coordinates ρ = ρ i ( z ) = φ i ( z ) and σ = σ j ( ζ ) = φ j ( ζ ) . (43)For notational clarity, ( ρ , σ ) denotes a point in the newcoordinate system, whereas ρ i ( z ) = φ i ( z ) and σ j ( ζ ) = φ j ( ζ ) are the respective transformations.Further, define ¯ K ( ρ i ( z ) , σ j ( ζ )) = λ j ( ζ ) K ( z, ζ ) . (44) EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 6
Then, differentiating (44) twice w. r. t. z and ζ as well asinserting the result in (41a) and (42a) yields λ i ( z )( φ (cid:48) i ( z )) ¯ K ρρ ( ρ , σ ) − λ j ( ζ )( φ (cid:48) j ( ζ )) ¯ K σσ ( ρ , σ )+ λ i ( z ) φ (cid:48)(cid:48) i ( z ) ¯ K ρ ( ρ , σ ) − λ j ( ζ ) φ (cid:48)(cid:48) j ( ζ ) ¯ K σ ( ρ , σ )= λ j ( ζ ) B [ K ]( z, ζ ) (45)after a simple calculation. Hence, one obtains the standardform ¯ K ρρ ( ρ , σ ) − ¯ K σσ ( ρ , σ ) + λ i ( z ) φ (cid:48)(cid:48) i ( z ) ¯ K ρ ( ρ , σ ) − λ j ( ζ ) φ (cid:48)(cid:48) j ( ζ ) ¯ K σ ( ρ , σ ) = λ j ( ζ ) B [ K ]( z, ζ ) (46)of the kernel PIDE, if λ i ( z )( φ (cid:48) i ( z )) = 1 , z ∈ (0 , , i = 1 , , . . . , n (47)holds. It is straightforward to verify that a solution of (47) is φ i ( z ) = (cid:90) z d ζ (cid:112) λ i ( ζ ) , z ∈ [0 , , i = 1 , , . . . , n. (48)Note that (48) also determines the function φ j ( ζ ) by appro-priate substitutions. The function φ i is invertible, because it isstrictly monotonically increasing. Hence, the inverse changeof coordinates z = z i ( ρ ) = φ − i ( ρ ) (49a) ζ = ζ j ( σ ) = φ − j ( σ ) (49b)related to (43) exists. Thereby, ( z, ζ ) denotes a point in theoriginal coordinate system, whereas z i ( ρ ) = φ − i ( ρ ) and ζ j ( σ ) = φ − j ( σ ) are the corresponding transformations.
2) Elimination of the first-order partial derivatives:
In viewof (48) the result φ (cid:48)(cid:48) i ( z ) = − λ (cid:48) i ( z )2 λ i ( z ) (cid:112) λ i ( z ) (50)can be deduced by differentiation. With this, the kernel PIDEs(46) can be rewritten as ¯ K ρρ ( ρ , σ ) − ¯ K σσ ( ρ , σ ) − κ i ( z i ( ρ )) ¯ K ρ ( ρ , σ ) + κ j ( ζ j ( σ )) ¯ K σ ( ρ , σ )= λ j ( ζ j ( σ )) B [ K ]( z i ( ρ ) , ζ j ( σ )) , (51)in which κ i ( z ) = λ (cid:48) i ( z )2 (cid:112) λ i ( z ) (52)and κ j ( ζ ) follows from the corresponding substitutions. Inview of (41a), (42a) and (51) the nonlinear transformation(43) introduces first-order partial derivatives in the kernelPIDEs. They, however, can readily be eliminated, because thecorresponding coefficient functions κ i ( z i ( ρ )) and κ j ( ζ j ( σ )) depend only on a single variable. Towards this end, the changeof coordinates ¯ K ( ρ , σ ) = ψ i ( ρ ) ψ j ( σ ) (cid:101) KKK ( ρ , σ ) (53)is introduced. From this and (44) the relation KKK ( z i ( ρ ) , ζ j ( σ )) = ψ i ( ρ ) ψ j ( σ ) λ j ( ζ j ( σ )) (cid:101) KKK ( ρ , σ ) (54) can be deduced. A straightforward but lengthy calculationshows that (51) takes the form (cid:101) KKK ρρ ( ρ , σ ) − (cid:101) KKK σσ ( ρ , σ ) = (cid:101) BBB [ (cid:101) K ]( ρ , σ ) (55)with (cid:101) BBB [ (cid:101) K ]( ρ , σ ) = (cid:16) a (cid:0) z i ( ρ ) , ζ j ( σ ) (cid:1) + µ c (cid:17) (cid:101) KKK ( ρ , σ ) (56) − λ j ( ζ j ( σ )) ψ i ( ρ ) ψ j ( σ ) F ( z i ( ρ ) , ζ j ( σ )) (cid:124) (cid:123)(cid:122) (cid:125) c ( ρ , σ ) + n (cid:88) k =1 (cid:101) K ik ( ρ , σ k ( ζ j ( σ ))) A kj ( ζ j ( σ )) λ j ( ζ j ( σ )) ψ k ( σ ) λ k ( ζ j ( σ )) ψ j ( σ ) (cid:124) (cid:123)(cid:122) (cid:125) c kj ( σ ) + (cid:90) z i ( ρ ) ζ j ( σ ) (cid:101) K ik ( ρ , σ k (¯ ζ )) F kj (¯ ζ, ζ j ( σ )) λ j ( ζ j ( σ )) ψ j ( σ ) ψ k ( σ k (¯ ζ )) λ k (¯ ζ ) (cid:124) (cid:123)(cid:122) (cid:125) c kj ( σ , ¯ ζ ) d¯ ζ and a ( z, ζ ) = − λ (cid:48)(cid:48) i ( z ) + λ (cid:48)(cid:48) j ( ζ ) + 3( λ (cid:48) i ( z )) λ i ( z ) − λ (cid:48) j ( ζ )) λ j ( ζ ) , (57)if ψ (cid:48) i ( ρ ) = κ i ( φ − i ( ρ ))2 ψ i ( ρ ) , ρ ∈ (0 , φ i (1)) , (58) i = 1 , , . . . , n , holds. A solution of (58) is ψ i ( ρ ) = (cid:115) λ i ( φ − i ( ρ )) λ i (0) (59)so that the transformation (53) is invertible. Note that (59) alsoyields ψ j ( σ ) by applying the corresponding substitutions.After utilizing the previous results, lengthy but straightfor-ward algebraic manipulations yield the BCs [ (cid:101) KKK ( ρ ,
0) = 0] λ i ≥ λ j , j ≤ m (60a) [ (cid:101) KKK σ ( ρ ,
0) + ( λ (cid:48) j (0)4 √ λ j (0) + q j (cid:112) λ j (0) (cid:124) (cid:123)(cid:122) (cid:125) c j ) (cid:101) KKK ( ρ , (cid:101) CCC [ (cid:101) K ]( ρ )] λ i ≥ λ j , j>m (60b)of the kernel PIDE (55) with (cid:101) CCC [ (cid:101) K ]( ρ ) = √ λ j (0) ψ i ( ρ ) A ( z i ( ρ )) (cid:124) (cid:123)(cid:122) (cid:125) c ( ρ ) (60c) − (cid:90) z i ( ρ )0 n (cid:88) k =1 (cid:101) K ik ( ρ , σ k (¯ ζ )) A ,kj (¯ ζ ) (cid:112) λ j (0) ψ k ( σ k (¯ ζ )) λ k (¯ ζ ) (cid:124) (cid:123)(cid:122) (cid:125) c kj (¯ ζ ) d¯ ζ and for i = j (cid:101) KKK ( ρ , ρ ) = − (cid:112) λ i (0) (cid:90) ρ (cid:0) A ( z i (¯ ρ )) + µ c (cid:1) d¯ ρ (cid:124) (cid:123)(cid:122) (cid:125) c i ( ρ ) . (61)Additionally, for i (cid:54) = j (cid:101) KKK (cid:0) ρ , σ j ( z i ( ρ )) (cid:1) = 0 (62a) EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 7 and (cid:101)
KKK ρ (cid:0) ρ , σ j ( z i ( ρ )) (cid:1) = (cid:113) λ i (0) λ j (0) λ i ( z i ( ρ )) λ j ( z i ( ρ )) λ j ( z i ( ρ )) − λ i ( z i ( ρ )) A ( z i ( ρ )) (cid:124) (cid:123)(cid:122) (cid:125) c ( ρ ) (62b)hold.
3) Canonical kernel equations:
The kernel PIDE (55) cannow be mapped into the canonical form . This only requires alinear change of coordinates so that no additional first-orderpartial derivatives are introduced in the result. The calculationof the related transformation coincides with the method forclassifying second-order PDEs (see, e. g., [17, Ch. 4]).The change of coordinates mapping the kernel PIDE (55)for (cid:101)
KKK ( ρ , σ ) = (cid:101) K ij ( ρ , σ ) into its canonical form depends onthe relation of the corresponding diffusion coefficients λ i and λ j . In particular, one obtains λ i ≥ λ j : ξ = ξ ρσij ( ρ , σ ) = ρ + σ (63a) η = η ρσij ( ρ , σ ) = ρ − σ (63b) λ i < λ j : ξ = ξ ρσij ( ρ , σ ) = φ i (1) + φ j (1) − ( ρ + σ ) (63c) η = η ρσij ( ρ , σ ) = − ( φ i (1) − φ j (1)) + ρ − σ (63d)with φ i defined in (48). In (63), ( ξ , η ) denotes a point in thecanonical coordinate system, whereas ξ ρσij ( ρ , σ ) and η ρσij ( ρ , σ ) are the corresponding transformations. Thereby, the superscript ρσ is introduced to identify the considered original coordinates.Besides mapping the kernel PIDE in its canonical form, thecoordinates (63c) and (63d) additionally lead to spatial domainsfor λ i < λ j , that coincide with the spatial domains in the case λ i ≥ λ j . This simplifies the subsequent derivation of the kernelintegral equations. By adding and subtracting the equations for ξ and η in (63), the inverse change of coordinates λ i ≥ λ j : ρ = ρ ij ( ξ , η ) = ξ + η (64a) σ = σ ij ( ξ , η ) = ξ − η (64b) λ i < λ j : ρ = ρ ij ( ξ , η ) = φ i (1) − ξ − η (64c) σ = σ ij ( ξ , η ) = φ j (1) − ξ + η (64d)is readily found. Note that a point ( ρ , σ ) can either becomputed by the functions ρ i ( z ) , σ j ( ζ ) (see (43)) from theoriginal coordinates, or by the functions ρ ij ( ξ , η ) and σ ij ( ξ , η ) according to (64) from the canonical coordinates. In order tosimplify the notation, introduce s = s ij = (cid:40) , λ i ≥ λ j − , λ i < λ j . (65)Insert (43) in (63) and use (65) to obtain the overall nonlinearcoordinate transformation ξ = ξ ij ( z, ζ ) = (1 − s )( φ i (1) + φ j (1)) + s ( φ i ( z ) + φ j ( ζ )) (66a) η = η ij ( z, ζ ) = − (1 − s )( φ i (1) − φ j (1)) + φ i ( z ) − φ j ( ζ ) (66b) from the original ( z, ζ ) -coordinates into the canonical coor-dinates. This transformation can be solved for ( z, ζ ) , whichyields z = z ij ( ξ , η ) = φ − i ( ( sξ + η ) + (1 − s ) φ i (1)) (67a) ζ = ζ ij ( ξ , η ) = φ − j ( ( sξ − η ) + (1 − s ) φ j (1)) . (67b)By combining the transformations (66) and (67), the canon-ical coordinates of different indices i and j can always betransformed into each other, which is necessary due to thecoupling between different kernel PIDEs. In particular, thetransformations ξ ξηkl ( ξ , η ) = ξ kl (cid:16) z ij ( ξ , η ) , ζ ij ( ξ , η ) (cid:17) (68a) η ξηkl ( ξ , η ) = η kl (cid:16) z ij ( ξ , η ) , ζ ij ( ξ , η ) (cid:17) (68b)can be introduced to transform from one canonical coordinatesystem ( ξ , η ) = ( ξ ij , η ij ) into the canonical coordinate systemof a different element ( ξ kl , η kl ) . Observe that the canonicalcoordinates can now be calculated from the original coordinatesby (66), from the ( ρ , σ ) -coordinates by (63) or from thecanonical coordinates of a different kernel element by (68). Anoverview of the various introduced coordinate transformationsis given in Figure 1.Now define G ( ξ , η ) = (cid:101) KKK ( ρ ij ( ξ , η ) , σ ij ( ξ , η )) , which canbe expressed as G ( ξ ρσij ( ρ , σ ) , η ρσij ( ρ , σ )) = (cid:101) KKK ( ρ , σ ) in thestandard form coordinates, differentiate it twice w. r. t. ρ and σ and insert it with (63) in (55). Using (65) provides thecompact notation of the canonical kernel PIDE s G ξη ( ξ , η ) = ˘ B [ G ]( ξ , η ) (69)after straightforward intermediate computations. Thereby, theabbreviation ˘ B [ G ]( ξ , η ) = (cid:16) a (cid:0) z ij ( ξ , η ) , ζ ij ( ξ , η ) (cid:1) + µ c (cid:17) G ( ξ , η ) − c ( ρ ij ( ξ , η ) , σ ij ( ξ , η ))+ n (cid:88) k =1 (cid:18) G ik ( ξ ξηik ( ξ , η ) , η ξηik ( ξ , η )) c kj ( σ ij ( ξ , η ))+ (cid:90) z ij ( ξ , η ) ζ ij ( ξ , η ) G ik (cid:16) ξ ik (cid:0) z ij ( ξ , η ) , ¯ ζ (cid:1) , η ik (cid:0) z ij ( ξ , η ) , ¯ ζ (cid:1)(cid:17) · c kj ( σ ij ( ξ , η ) , ¯ ζ )d¯ ζ (cid:19) (70)was used (see (56)). Remark 5:
It should be noted that the proposed three steptransformation into the kernel PIDE (69) is much simpleras directly applying the classical change of coordinates forclassifying second-order PDEs to the kernel PIDEs (41a) and(42a). This is due to the fact that the latter approach yieldsfirst-order partial derivatives with coefficient functions whichexhibit an involved structure. Consequently, the determinationof the corresponding transformation for their elimination isimpeded. (cid:47)
Remark 6:
For mutually different constant diffusion coeffi-cients λ i = const. , i = 1 , , . . . , n , the overall linear change EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 8 ( z , ζ ) ( ρ , σ ) ( ξ , η ) ρ i ( z ) , σ j ( ζ ) z i ( ρ ) , ζ j ( σ ) ξ ρσij ( ρ , σ ) , η ρσij ( ρ , σ ) ρ ij ( ξ , η ) , σ ij ( ξ , η ) ξ ij ( z, ζ ) , η ij ( z, ζ ) z ij ( ξ , η ) , ζ ij ( ξ , η ) ij → klξ ξηkl ( ξ , η ) η ξηkl ( ξ , η ) original coordinates ( z, ζ ) standard form coordinates ( ρ , σ ) canonical coordinates ( ξ , η ) Fig. 1. Overview of the changes of coordinates for mapping the kernel equations into their canonical form. of coordinates λ i ≥ λ j : ξ = z √ λ i + ζ √ λ j and η = z √ λ i − ζ √ λ j (71a) λ i < λ j : ξ = √ λ i + √ λ j − ( z √ λ i + ζ √ λ j ) (71b) η = − ( √ λ i + √ λ j ) + z √ λ i − ζ √ λ j (71c)can be applied to obtain the canonical kernel PIDE in a singlestep. This directly follows from (43), (48) and (63). (cid:47) For the derivation of the canonical kernel equations itremains to determine the BCs for (69). They result fromapplying (63) and (64) to (60)–(62). For this, the spatial domainof (69) is investigated.Consider the boundary ζ = z , z ∈ [0 , , of the spatialdomain w. r. t. the original kernel PIDEs (41a) and (42a) andsubstitute ζ = z in (66) to obtain ξ = (1 − s )( φ i (1) + φ j (1)) + sβ ( z ) (72a) η = − (1 − s )( φ i (1) − φ j (1)) + φ i ( z ) − φ j ( z ) (72b)with β ( z ) = φ i ( z ) + φ j ( z ) . The function β ( z ) is strictlymonotonically increasing, because φ i and φ j have this property(see (48)). Hence, the inverse β − exists and thus (72a) canbe solved for z , resulting in z l ( ξ ) = β − (cid:16) sξ + (1 − s )( φ i (1) + φ j (1)) (cid:17) . (73)Inserting this into (72b) yields the lower boundary η l ( ξ ) = − (1 − s )( φ i (1) − φ j (1)) + φ i ( z l ( ξ )) − φ j ( z l ( ξ )) (74)of the domain D ij in Figure 2. It is easy to show that for i (cid:54) = j , η (cid:48) l ( ξ ) < for ξ ∈ [0 , φ i (1) + φ j (1)] so that η l ( ξ ) is strictly monotonically decreasing. This implies that theresulting boundary of (69) only evolves in the fourth quadrantof the ( ξ , η ) -plane (see Figure 2). In contrast, inserting i = j in (74) shows that [ η l ( ξ ) = 0] i = j (75)holds so that the spatial domain is the triangle η ∈ (0 , φ i (1)) , ξ ∈ ( η , φ i (1) − η ) . By a similar but very elementary reasoningone can transform the remaining boundaries of the kernel ξη a φ i (1) + φ j (1) ba D ij Γ Γ : η l ( ξ ) Fig. 2. Spatial domains D ij of the canonical kernel equations (76) with λ i ≥ λ j : a = φ i (1) , b = φ i (1) − φ j (1) and λ i < λ j : a = φ j (1) , b = − ( φ i (1) − φ j (1)) . The bold lines at Γ and Γ represent the BCs of thekernel BVP, in which Γ is an artificial BC for λ i < λ j (see (76e)). Thedash dotted arrows (’ − · − ’) indicate the integration in the ξ -direction andthe solid arrow (’–’) the integration in the η -direction. Note that in the case i = j the spatial domain is the triangle η ∈ (0 , φ i (1)) , ξ ∈ ( η , φ i (1) − η ) ,as b = 0 and η l ( ξ ) ≡ . PIDEs (41a) and (42a) into the ( ξ , η ) -coordinates with (66).The resulting spatial domain D ij is depicted in Figure 2.With the spatial domain of the canonical coordinates known,the BCs (60)–(62) can be mapped into ( ξ , η ) -coordinates bymaking use of (63) and (64), resulting in the canonical kernelequations s G ξη ( ξ , η ) = ˘ B [ G ]( ξ , η ) (76a) [ G ( η , η ) = 0] λ i ≥ λ j , j ≤ m (76b) [ G ξ ( η , η ) − G η ( η , η ) + c j G ( η , η ) = ˘ C [ G ]( η )] λ i ≥ λ j j >m (76c) [ G ( ξ ,
0) = c i ( ξ )] i = j (76d) [ G ( η , η ) = g f ( η )] λ i <λ j (76e) [ G ( ξ , η l ( ξ )) = 0] i (cid:54) = j (76f) [ sG ξ ( ξ , η l ( ξ )) + G η ( ξ , η l ( ξ )) = c (cid:0) ρ ij ( ξ , η l ( ξ )) (cid:1) ] i (cid:54) = j (76g)for i, j = 1 , , . . . , n , in which (76a) is defined on D ij (see EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 9
Figure 2) and ˘ C [ G ]( η ) = c ( ρ ij ( η , η )) (77) − z ij ( η , η ) (cid:90) n (cid:88) k =1 G ik (cid:16) ξ ik (cid:0) z ij ( η , η ) , ¯ ζ (cid:1) , η ik (cid:0) z ij ( η , η ) , ¯ ζ (cid:1)(cid:17) c kj (¯ ζ )d¯ ζ was used (see (60c)). The BC (76e) at Γ is an artificial BC ,that has to be introduced in order to ensure well-posed kernelequations. Thereby, g f ∈ C [0 , φ j (1)] is a degree of freedomand can be chosen arbitrarily. Remark 7:
It is noteworthy that two BCs are needed at theboundary Γ , because the derivation of the integral equationsrequires one BC at Γ for the integration in the direction of ξ and one BC at Γ for the integration in the direction of η (seeFigure 2). Furthermore, the BCs at Γ and Γ utilized for theintegration w. r. t. the ξ -direction imply that the ξ -axis within D ij constitutes a line of separation defining two differentparts of the corresponding kernel on the spatial domain D ij ,leading to piecewise defined kernels for i (cid:54) = j . It can be shownthat the line of separation is no longer a straight line in theoriginal coordinates but a strictly monotonically increasingcurve defined by ζ = φ − j ( φ i ( z )) , z ∈ [0 , , for λ i > λ j and z = φ − i ( φ i (1) − φ j (1) + φ j ( ζ )) , ζ ∈ [0 , , if λ i < λ j . (cid:47) C. Kernel Integral Equations
In order to solve the canonical kernel equations (76) utilizinga successive approximation, they are transformed into integralequations. The convergence analysis can be simplified byrewriting the second-order PIDE (76a) into a system of twofirst-order PIDEs. By introducing H ( ξ , η ) = G ξ ( ξ , η ) , (78)the PIDE (76a) can be written as G ξ ( ξ , η ) = H ( ξ , η ) (79a) H η ( ξ , η ) = s ˘ BBB [ G ]( ξξξ, ηηη ) . (79b)To formulate the corresponding BCs, solve (76c) for G η ( η , η ) ,insert the result into d η G ( η , η ) = G ξ ( η , η ) + G η ( η , η ) (80)and use (78). Then, an integration w. r. t. η yields (cid:104) G ( η , η ) = η (cid:90) (cid:16) H (¯ η, ¯ η ) + c j G (¯ η, ¯ η ) − ˘ CCC [ G ](¯ η ) (cid:17) d¯ η (cid:105) λ i ≥ λ j j >m . (81)In view of (76f) and (78) one obtains d ξ G ( ξ , η l ( ξ )) = H ( ξ , η l ( ξ )) + G η ( ξ , η l ( ξ )) η (cid:48) l ( ξ ) = 0 . (82)Solving (76g) for G η ( ξ , η l ( ξ )) and inserting the result in (82)yields (cid:104) H ( ξ , η l ( ξ )) = ccc (cid:0) ρ ij ( ξ , η l ( ξ )) (cid:1) η (cid:48) l ( ξ ) sη (cid:48) l ( ξ ) − (cid:124) (cid:123)(cid:122) (cid:125) ccc ( ξ ) (cid:105) i (cid:54) = j (83) after a simple rearrangement. It can easily be shown that − < η (cid:48) l ( ξ ) < for ξ ∈ [0 , φ i (1) + φ j (1)] . Hence, the denominatorin (83) cannot be zero. Finally, differentiate (76d) w. r. t. ξ anduse (78), the definition of c i in (61) and (75) to obtain (cid:2) H ( ξ ,
0) = H ( ξ , η l ( ξ )) = − √ λ i (0)4 (cid:0) AAA ( z ij ( ξ , µ c (cid:1)(cid:124) (cid:123)(cid:122) (cid:125) ccc ( ξ ) (cid:3) i = j . (84)For the ease of presentation, the left boundary of the spatialdomain D ij can be introduced as ξ l ( η ) = (cid:40) η , η ≥ η − l ( η ) , η < (85)with η − l being the inverse of (74). With this, an alternativeway to write (76f) is [ GGG ( ξ l ( η ) , η ) = 0] i (cid:54) = j, η < . (86)Now, the BC on the left boundary of the spatial domain D ij can be compactly written as G ( ξ l ( η ) , η ) = (cid:40) G ( η , η ) , η ≥ , i (cid:54) = j, η < (87)in view of (85) and (86). With this and G ( η , η ) determinedby (76b), (76e) and (81), the BCs needed for the derivation ofthe integral equations are G ( ξ l ( η ) , η ) = (cid:82) η (cid:0) H (¯ η, ¯ η ) + c j G (¯ η, ¯ η ) − ˘ CCC [ G ](¯ η ) (cid:1) d¯ η,λ i ≥ λ j , j > m, η ≥ , λ i ≥ λ j , j ≤ m, η ≥ ggg f ( η ) , λ i < λ j , η ≥ , i (cid:54) = j, η < (88a)and H ( ξ , η l ( ξ )) = (cid:40) ccc ( ξ ) , i = jccc ( ξ ) , i (cid:54) = j . (88b)Now, (79a) is integrated w. r. t. ξ and (79b) w. r. t. η , leadingto G ( ξ , η ) = G ( ξ l ( η ) , η ) + (cid:90) ξξ l ( η ) H ( ¯ ξ, η )d ¯ ξ (89a) H ( ξ , η ) = H ( ξ , η l ( ξ )) + (cid:90) ηη l ( ξ ) 14 s ˘ B [ G ]( ξξξ, ¯ η )d¯ η, (89b)when taking the boundaries of D ij into account. After inserting(88a) with (77) into (89a), the first kernel integral equation reads G ( ξ , η ) = G ( η ) + FFF G [ G, H ]( ξ , η ) = (cid:101) FFF G [ G, H ]( ξ , η ) , (90)where G ( η ) = − (cid:82) η ccc ( ρ ij (¯ η, ¯ η ))d¯ η, λ i ≥ λ j , j > m, η ≥ ggg f ( η ) , λ i < λ j , η ≥ , else (91) EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 10 and
FFF G [ G, H ]( ξ , η ) = (cid:90) ξξ l ( η ) H ( ¯ ξ, η )d ¯ ξ (92) + (cid:104) (cid:90) η (cid:16) H (¯ η, ¯ η ) + c j G (¯ η, ¯ η ) + (cid:90) z ij (¯ η, ¯ η )0 n (cid:88) k =1 c kj (¯ ζ ) · G ik (cid:0) ξ ik ( z ij (¯ η, ¯ η ) , ¯ ζ ) , η ik ( z ij (¯ η, ¯ η ) , ¯ ζ ) (cid:1) d¯ ζ (cid:17) d¯ η (cid:105) λ i ≥ λ j j >m η > hold.Furthermore, after inserting (70) in (89b) the second kernelintegral equation follows as H ( ξ , η ) = H ( ξ ) + FFF H [ G ]( ξ , η ) = (cid:101) FFF H [ G ]( ξ , η ) (93)with H ( ξ ) = H ( ξ , η l ( ξ )) − s (cid:90) ηη l ( ξ ) ccc ( ρ ij ( ξ , ¯ η ) , σ ij ( ξ , ¯ η ))d¯ η (94)and H ( ξ , η l ( ξ )) according to (88b) as well as FFF H [ G ]( ξ , η ) =14 s (cid:90) ηη l ( ξ ) (cid:18)(cid:0) a (cid:0) z ij ( ξ , ¯ η ) , ζ ij ( ξ , ¯ η ) (cid:1) + µ c (cid:1) G ( ξ , ¯ η )+ n (cid:88) k =1 (cid:16) G ik (cid:0) ξ ξηik ( ξ , ¯ η ) , η ξηik ( ξ , ¯ η ) (cid:1) c kj ( σ ij ( ξ , ¯ η ))+ (cid:90) z ij ( ξ , ¯ η ) ζ ij ( ξ , ¯ η ) c kj (cid:0) σ ij ( ξ , ¯ η ) , ¯ ζ (cid:1) · G ik (cid:0) ξ ik (cid:0) z ij ( ξ , ¯ η ) , ¯ ζ (cid:1) , η ik (cid:0) z ij ( ξ , ¯ η ) , ¯ ζ (cid:1)(cid:1) d¯ ζ (cid:17)(cid:19) d¯ η. (95) D. Successive Approximation
To compute the solution of the kernel integral equations (90)and (93), the method of successive approximations presented in[21] is extended to the considered multivariable case. However,there is not only a single integral equation, but n integralequations for each G and H , which are coupled.With the recursions G l +1 = (cid:101) FFF G [ G l , H l ] , l ∈ N (96a) H l +1 = (cid:101) FFF H [ G l ] , G = 0 , H = 0 , (96b)implied by (90) and (93), the corresponding fixpoints G = lim N →∞ G N and H = lim N →∞ H N (97)are the solutions of the integral equations, if they exist. Thelatter amounts to proving the corresponding convergence, forwhich the equivalent representation G = ∞ (cid:88) l =0 ∆ G l and H = ∞ (cid:88) l =0 ∆ H l (98)with ∆ G l +1 = FFF G [∆ G l , ∆ H l ] , ∆ G = G ( η ) (99a) ∆ H l +1 = FFF H [∆ G l ] , ∆ H = HHH ( ξ ) (99b) is considered.As the integral operators (92) and (95) contain sums overthe kernel elements that are defined in different coordinatesystems, a growth assumption is needed, which is independentfrom the coordinate systems of the kernel elements. For thispurpose, assume | ∆ G l ( ξ ij ( z, ζ ) , η ij ( z, ζ )) | ≤ M l +1 l ! ( z − γζ ) l (100a) | ∆ H l ( ξ ij ( z, ζ ) , η ij ( z, ζ )) | ≤ M l +1 l ! ( z − γζ ) l (100b)with γ ∈ (cid:32) max λ i <λ j (cid:115) λ i ( zzz ∆ ) λ j ( zzz ∆ ) , (cid:33) . (101)Thereby, zzz ∆ = argmin( | λ i ( z ) − λ j ( z ) | ) (102)is the point of minimal difference between the diffusioncoefficients λ i and λ j .The growth assumptions (100) can be inserted in (98) andexpressed in the canonical coordinates with (67). This gives | G ( ξ , η ) | ≤ ∞ (cid:88) l =0 M l +1 l ! ( z ij ( ξ , η ) − γζ ij ( ξ , η )) l (103a) = M e M ( z ij ( ξ , η ) − γζ ij ( ξ , η )) = M e M ( z − γζ ) | H ( ξ , η ) | ≤ ∞ (cid:88) l =0 M l +1 l ! ( z ij ( ξ , η ) − γζ ij ( ξ , η )) l (103b) = M e M ( z ij ( ξ , η ) − γζ ij ( ξ , η )) = M e M ( z − γζ ) showing that the series (98) converge absolutely and uniformly.Hence, the piecewise continuous solution of the kernel integralequations (90) and (93) can be computed via (96), if (100) isverified for all elements ∆ G l and ∆ H l .Considering the initial values ∆ G and ∆ H given by (91)and (94), it can easily be seen that (100) is valid for l = 0 , ifthe coefficient functions ccc , ccc , ccc and ccc , as well as the degreeof freedom ggg f are bounded. The definitions of the coefficientfunctions directly show that this is always true, if the systemparameters Λ , A , A and F in (1) are bounded. What remainsto be proven is that the integral operators (92) and (95) usedfor the update law (99) preserve the growth assumption (100)for all elements ∆ G l and ∆ H l , l > .
1) Preliminaries:
To prove the validity of the growthassumptions (100) for l > , it is convenient to take somepreliminary algebraic considerations first.The inverse transformation (67) can be inserted in z − γζ and differentiated w. r. t. the canonical coordinates ( ξ , η ) . Toshorten the presentation, the points ( z , ζ ) will be writteninstead of the functions z ij ( ξ , η ) and ζ ij ( ξ , η ) in the following.This yields ∂ ξ ( z − γ ζ ) = ∂ ξ (cid:0) φ − i ( ( sξ + η ) + 12 (1 − s ) φ i (1)) (cid:1) (104) − γ∂ ξ (cid:0) φ − j ( ( sξ − η ) + (1 − s ) φ j (1)) (cid:1) . EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 11
Considering the well-known result for the derivative of theinverse function of (48) dd y φ − i ( y ) (cid:12)(cid:12) y = φ i ( z ) = 1 φ (cid:48) i ( z ) = (cid:112) λ i ( z ) , i = 1 , , . . . , n (105)and inserting (67) leads to the result ∂ ξ ( z − γ ζ ) = s (cid:16)(cid:113) λ i ( z ) − γ (cid:113) λ j ( ζ ) (cid:17) . (106)With the definitions (65) and (101) it is easy to see that (106)is strictly positive. Repeating this procedure for the derivativew. r. t. η shows that in summary ∂ ξ ( z − γ ζ ) > (107a) ∂ η ( z − γ ζ ) > (107b)hold. In view of the spatial domain ≤ ζ ≤ z ≤ in theoriginal coordinates and (101), the inequality ≤ z − γζ ≤ (108)can be deduced. With this and (107), the relations ≤ z ij ( ξ ∗ , η ) − γζ ij ( ξ ∗ , η ) z ij ( ξ , η ) − γζ ij ( ξ , η ) ≤ (109a) ≤ z ij ( ξ , η ∗ ) − γζ ij ( ξ , η ∗ ) z ij ( ξ , η ) − γζ ij ( ξ , η ) ≤ (109b)are satisfied for all points ξ ∗ ≤ ξ and η ∗ ≤ η .Moreover, calculating the derivative of ( z − γ ζ ) l +1 w. r. t. ξ and using (106) results in ∂ ξ ( z − γ ζ ) l +1 = ( l + 1)( z − γ ζ ) l ∂ ξ ( z − γ ζ ) (110) = ( l + 1)( z − γ ζ ) l s (cid:16)(cid:113) λ i ( z ) − γ (cid:113) λ j ( ζ ) (cid:17) ≥ ( l + 1)( z − γ ζ ) l s (cid:16)(cid:113) λ i ( z ∆ ) − γ (cid:113) λ j ( z ∆ ) (cid:17)(cid:124) (cid:123)(cid:122) (cid:125)(cid:12)(cid:12) √ λ i ( z ∆ ) − γ √ λ j ( z ∆ ) (cid:12)(cid:12) . At this point, the importance of the choice of γ according to(101) becomes visible. It ensures that the right sides of (106)and (110) are positive which is the basis for (109) and thefollowing bound estimations.Solving (110) for ( z − γ ζ ) l and integrating both sides w. r. t. ξ leads to the upper bound (cid:90) ξξ l (cid:0) z ij ( ¯ ξ, η ) − γζ ij ( ¯ ξ, η ) (cid:1) l d ¯ ξ ≤ (cid:104)(cid:0) z ij ( ¯ ξ, η ) − γζ ij ( ¯ ξ, η ) (cid:1) l +1 (cid:105) ξξ l ( l + 1) (cid:12)(cid:12)(cid:112) λ i ( z ∆ ) − γ (cid:112) λ j ( z ∆ ) (cid:12)(cid:12) . (111a)With an identical reasoning, it is possible to prove the inequal-ity (cid:90) ηη l (cid:0) z ij ( ξ , ¯ η ) − γζ ij ( ξ , ¯ η ) (cid:1) l d¯ η ≤ (cid:104)(cid:0) z ij ( ξ , ¯ η ) − γζ ij ( ξ , ¯ η ) (cid:1) l +1 (cid:105) ηη l ( l + 1) (cid:0)(cid:112) λ i ( z Σ ) + γ (cid:112) λ j ( z Σ ) (cid:1) (111b) and similarly for λ i ≥ λ j (cid:90) η (cid:0) z ij (¯ η, ¯ η ) − γζ ij (¯ η, ¯ η ) (cid:1) l d¯ η ≤ (cid:104)(cid:0) z ij (¯ η, ¯ η ) − γζ ij (¯ η, ¯ η ) (cid:1) l +1 (cid:105) η ( l + 1) (cid:112) λ i ( z Σ ) , (111c)where z Σ = argmin( λ i ( z ) + λ j ( z )) (112)is the point for which the sum of the diffusion coefficients λ i and λ j is minimal.By making use of (66) and (67), it can be verified that z ik (cid:0) ξ ik ( z ij (¯ η, ¯ η ) , ¯ ζ ) , η ik ( z ij (¯ η, ¯ η ) , ¯ ζ ) (cid:1) = z ij (¯ η, ¯ η ) (113a) ζ ik (cid:0) ξ ik ( z ij (¯ η, ¯ η ) , ¯ ζ ) , η ik ( z ij (¯ η, ¯ η ) , ¯ ζ ) (cid:1) = ¯ ζ (113b)is satisfied, which will significantly simplify the evaluation ofthe sums in the following steps.
2) Proof of the growth assumptions for l > : Now thegrowth assumptions (100) can be inserted into the fixpointiteration (99) with the operators (92) and (95), to prove itsvalidity via induction.With (100), (113), the triangle inequality and factoring M l +1 l ! out, the recursion (99a) and (92) give | ∆ G l +1 ( ξ , η ) | ≤ M l +1 l ! (cid:16) (cid:90) ξξ l (cid:16) z ij ( ¯ ξ, η ) − γζ ij ( ¯ ξ, η ) (cid:17) l d ¯ ξ + (cid:104) (cid:90) η (cid:16) (2 + c j ) (cid:0) z ij (¯ η, ¯ η ) − γζ ij (¯ η, ¯ η ) (cid:1) l (114) + (cid:90) z ij (¯ η, ¯ η )0 n (cid:88) k =1 (cid:0) z ij (¯ η, ¯ η ) − γ ¯ ζ (cid:1) l c kj (¯ ζ )d¯ ζ (cid:17) d¯ η (cid:105) λ i ≥ λ j j >m η > (cid:17) . Due to the fact that the growth assumptions (100) wereformulated in the original coordinates, which led to thesimplification (113), the term ( z ij (¯ η, ¯ η ) − γ ¯ ζ ) l can now befactored out of the sum, i. e., n (cid:88) k =1 ( z ij (¯ η, ¯ η ) − γ ¯ ζ ) l c kj (¯ ζ ) = ( z ij (¯ η, ¯ η ) − γ ¯ ζ ) l n (cid:88) k =1 c kj (¯ ζ ) ≤ ( z ij (¯ η, ¯ η ) − γ ¯ ζ ) l n ¯ c (115)is obtained, in which c kj (¯ ζ ) ≤ ¯ c , k, j = 1 , , . . . , n , ¯ ζ ∈ [0 , ,follows from the boundedness of the system parameters.With (115), the remaining inner integral in (114) can beevaluated as (cid:90) z ij (¯ η, ¯ η )0 (cid:16) z ij (¯ η, ¯ η ) − γ ¯ ζ (cid:17) l d¯ ζ (116) = − γ ( l + 1) (cid:104)(cid:16) z ij (¯ η, ¯ η ) − γ ¯ ζ (cid:17) l +1 (cid:105) z ij (¯ η, ¯ η )0 = 1 γ ( l + 1) (cid:16) z ij (¯ η, ¯ η ) l +1 − (cid:0) z ij (¯ η, ¯ η ) − γz ij (¯ η, ¯ η ) (cid:1) l +1 (cid:17) = 1 γ ( l + 1) (cid:16) z ij (¯ η, ¯ η ) − γζ ij (¯ η, ¯ η ) (cid:17) l +1 · (cid:32)(cid:18) z ij (¯ η, ¯ η ) z ij (¯ η, ¯ η ) − γζ ij (¯ η, ¯ η ) (cid:19) l +1 − (cid:18) z ij (¯ η, ¯ η ) − γz ij (¯ η, ¯ η ) z ij (¯ η, ¯ η ) − γζ ij (¯ η, ¯ η ) (cid:19) l +1 (cid:33) . EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 12
As the square bracket in (114) only exists for λ i ≥ λ j , theinverse transformation (67) has to be considered for s = 1 in(116) (see (65)). Consequently, inserting ξ = η = ¯ η in (67b)results in ζ ij (¯ η, ¯ η ) = 0 . Thus, (116) can be simplified to (cid:90) z ij (¯ η, ¯ η )0 (cid:16) z ij (¯ η, ¯ η ) − γ ¯ ζ (cid:17) l d¯ ζ ≤ γ ( l + 1) (cid:0) z ij (¯ η, ¯ η ) − γζ ij (¯ η, ¯ η ) (cid:1) l +1 (cid:0) − (1 − γ ) l +1 (cid:1) ≤ γ ( l + 1) (cid:0) z ij (¯ η, ¯ η ) − γζ ij (¯ η, ¯ η ) (cid:1) l +1 ≤ γ (cid:0) z ij (¯ η, ¯ η ) − γζ ij (¯ η, ¯ η ) (cid:1) l , (117)because of (101), (108) and l > . The first appearance of ζ ij (¯ η, ¯ η ) is kept deliberately, because it removes the necessityof case distinctions in the sequel. Inserting (117) in (114) yields | ∆ G l +1 ( ξ , η ) | ≤ M l +1 l ! (cid:16) (cid:90) ξξ l (cid:0) z ij ( ¯ ξ, η ) − γζ ij ( ¯ ξ, η ) (cid:1) l d ¯ ξ + (cid:104) (2 + c j + n ¯ cγ ) (cid:124) (cid:123)(cid:122) (cid:125) ¯¯ c (cid:90) η (cid:0) z ij (¯ η, ¯ η ) − γζ ij (¯ η, ¯ η ) (cid:1) l d¯ η (cid:105) λ i ≥ λ j j >m η > (cid:17) . (118)To evaluate the integrals in (118), insert (111a) and (111c).After factoring out l +1 ( z ij ( ξ , η ) − γζ ij ( ξ , η )) l +1 , the relation | ∆ G l +1 ( ξ , η ) | ≤ M l +1 l !( l + 1) (cid:0) z ij ( ξ , η ) − γζ ij ( ξ , η ) (cid:1) l +1 (cid:18) | (cid:112) λ i ( z ∆ ) − γ (cid:112) λ j ( z ∆ ) | (cid:16) − (cid:16) z ij ( ξ l , η ) − γζ ij ( ξ l , η ) z ij ( ξ , η ) − γζ ij ( ξ , η ) (cid:17) l +1 (cid:17) + (cid:104) ¯¯ c (cid:112) λ i ( z Σ ) (cid:16)(cid:16) z ij ( η , η ) − γζ ij ( η , η ) z ij ( ξ , η ) − γζ ij ( ξ , η ) (cid:17) l +1 − (cid:16) z ij (0 , − γζ ij (0 , z ij ( ξ , η ) − γζ ij ( ξ , η ) (cid:124) (cid:123)(cid:122) (cid:125) =0 (cid:17) l +1 (cid:17)(cid:105) λ i ≥ λ j j >m η > (cid:19) (119)is obtained. In order to prove convergence, it is necessary thatthe absolute values of all terms with the exponent l + 1 arelower or equal to one. Observe that η ≤ ξ in the domain D ij (see Figure 2) and obviously ξ l ≤ ξ as it is the left boundaryof the domain. Thus, (109a) holds for the fractions in (119),why it can be further simplified to | ∆ G l +1 ( ξ , η ) | ≤ M l +1 ( l + 1)! (cid:0) z ij ( ξ , η ) − γζ ij ( ξ , η ) (cid:1) l +1 · (cid:18) (cid:112) λ i ( z ∆ ) − γ (cid:112) λ j ( z ∆ ) + (cid:104) ¯¯ c (cid:112) λ i ( z Σ ) (cid:105) λ i ≥ λ j j >m η > (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) ≤ ˜ c ≤ M l +1 ˜ c ( l + 1)! (cid:0) z ij ( ξ , η ) − γζ ij ( ξ , η ) (cid:1) l +1 , (120)where ˜ c is a sufficiently large but finite number. Now an M ≥ ˜ c (121)can always be found such that | ∆ G l +1 ( ξ , η ) | ≤ M l +2 ( l + 1)! (cid:0) z ij ( ξ , η ) − γζ ij ( ξ , η ) (cid:1) l +1 (122) is satisfied. Comparing (122) with (100a) shows the validityof the growth assumption (100a) for all l > by induction.Now (100a) is inserted into the second update law (99b) withthe integral operator (95). After applying the simplification(113), factoring out M l +1 l ! and introducing the upper bound ˘ c of all appearing coefficients, the result | ∆ H l +1 ( ξ , η ) | ≤ M l +1 ˘ cl ! (cid:90) ηη l (cid:18)(cid:0) z ij ( ξ , ¯ η ) − γζ ij ( ξ , ¯ η ) (cid:1) l + n (cid:88) k =1 (cid:16)(cid:0) z ij ( ξ , ¯ η ) − γζ ij ( ξ , ¯ η ) (cid:1) l + (cid:90) z ij ( ξ , ¯ η ) ζ ij ( ξ , ¯ η ) (cid:0) z ij ( ξ , ¯ η ) − γ ¯ ζ (cid:1) l d¯ ζ (cid:17)(cid:19) d¯ η (123)follows. As all terms in the sum are independent of k , the sumcan be replaced by a multiplication with n . Moreover, the innerintegral can be evaluated in a similar way as in (116), whichshows that it can be replaced by γ ( z ij ( ξ , ¯ η ) − γζ ij ( ξ , ¯ η )) l in(123). In summary, the estimate | ∆ H l +1 ( ξ , η ) | (124) ≤ M l +1 ˘ cl ! η (cid:90) η l (2+ n + nγ ) (cid:0) z ij ( ξ , ¯ η ) − γζ ij ( ξ , ¯ η ) (cid:1) l d¯ η is obtained. For the evaluation of the integral, the relation(111b) is applied. Then ( z ij ( ξ , η ) − γζ ij ( ξ , η )) l +1 is factoredout after inserting the limits of the integral, which leads to | ∆ H l +1 ( ξ , η ) | ≤ M l +1 ˘ cl !( l + 1) (cid:0) z ij ( ξ , η ) − γζ ij ( ξ , η ) (cid:1) l +1 (125) n + nγ ) (cid:112) λ i ( z Σ ) + γ (cid:112) λ j ( z Σ ) (cid:16) − (cid:16) z ij ( ξ , η l ) − γζ ij ( ξ , η l ) z ij ( ξ , η ) − γζ ij ( ξ , η ) (cid:17) l +1 (cid:17) . It is clear that η l ≤ η , as it is the lower boundary of thedomain D ij . Thus, the relation (109b) can be used for furthersimplifying (125) to | ∆ H l +1 ( ξ , η ) | ≤ M l +2 ( l + 1)! (cid:0) z ij ( ξ , η ) − γζ ij ( ξ , η ) (cid:1) l +1 (126)for some M ≥ ˘ c n + nγ ) (cid:112) λ i ( z Σ ) + γ (cid:112) λ j ( z Σ ) . (127)Observe that with the definition of γ in (101) and therequirements for M , given by (121) and (127), an appropriate,bounded M can always be found. With (126) the validity of thegrowth assumption (100b) for all l > follows by induction.As a result, the kernel integral equations (90) and (93) admita piecewise continuous solution. The next theorem states thatthe corresponding kernel K ( z, ζ ) is also a piecewise classicalsolution of (24). Theorem 2 (Kernel equations):
The kernel equations (24)have a piecewise C -solution on the spatial domain ≤ ζ ≤ z ≤ . Proof:
As a result of the previous paragraph, the kernelintegral equations (90) and (93) admit a piecewise continuoussolution. The regularity of the kernel obtained from thesuccessive approximation remains to be verified. For this
EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 13 purpose, differentiate (90) twice w. r. t. η as well as (93) w. r. t. ξ and η . By utilizing the piecewise continuity of G ( ξ , η ) and H ( ξ , η ) as well as the assumed regularity of Λ , A and F , thisyields piecewise continuous functions. By direct substitution ofthe kernel integral equations in the respective kernel equationsit follows that (76) and consequently (24) have a piecewiseclassical solution, which proves the theorem.V. C LOSED - LOOP STABILITY
In order to prove the stability of the closed-loop systemresulting from applying (9) to (1), the bounded invertibility ofthe backstepping transformation (10) has to be verified. Forthis, the inverse backstepping transformation x ( z, t ) = ˜ x ( z, t ) + (cid:90) z L ( z, ζ )˜ x ( ζ, t )d ζ (128)with the integral kernel L ( z, ζ ) ∈ R n × n is introduced. Inthe proof of the next theorem it is demonstrated that thetransformation (128) mapping the target system (11) back intothe original coordinates exists. Thereby, the integral kernel L ( z, ζ ) is a piecewise C -function implying the boundednessof (128). Hence, the exponential stability of (11) implies thesame in the original coordinates, which is the result of thefollowing theorem. Theorem 3 (Closed-loop stability):
Assume that µ c > µ max (see Theorem 1). Then, the closed-loop system (1) and (9) is ex-ponentially stable in the L -norm (cid:107) h (cid:107) = ( (cid:82) (cid:107) h ( z ) (cid:107) C n d z ) / ,i. e., (cid:107) x ( t ) (cid:107) ≤ M e ( µ max − µ c ) t (cid:107) x (0) (cid:107) , t ≥ (129)for all x (0) ∈ ( H (0 , n compatible with the BCs of theclosed-loop system (1), (9) and an M ≥ . Proof:
The transformation (128) maps (11) into theoriginal coordinates, if L ( z, ζ ) is the solution of the inversekernel equations Λ( z ) L zz ( z, ζ ) − ( L ( z, ζ )Λ( ζ )) ζζ = E [ L ]( z, ζ ) (130a) L ( z, E = −F [ L ]( z ) E (130b) L ζ ( z, E + L ( z, (cid:48) (0) E + Λ(0) E Q ) (130c) = (cid:0) A ( z ) + F [ L ]( z ) (cid:1) E Λ( z ) L (cid:48) ( z, z )+Λ( z ) L z ( z, z )+ L ζ ( z, z )Λ( z ) + L ( z, z )Λ (cid:48) ( z )= − ( A ( z )+ µ c I ) (130d) L ( z, z )Λ( z ) − Λ( z ) L ( z, z ) = 0 (130e) L (0 ,
0) = 0 (130f)with (130a) defined on < ζ < z < and E [ L ]( z, ζ ) = − ( A ( z ) + µ c I ) L ( z, ζ ) − F ( z, ζ ) − (cid:90) zζ F ( z, ζ (cid:48) ) L ( ζ (cid:48) , ζ )d ζ (cid:48) (131a) F [ L ]( z ) = (cid:101) A ( z ) + (cid:90) z L ( z, ζ ) (cid:101) A ( ζ )d ζ. (131b)They follow from the same reasoning as in Section III-B. Note,that therein (cid:101) A ( z ) depends on the kernel K ( z, ζ ) (see (31)) andthus is assumed to be known. Considering the component formof (130), the BCs originating from (130b) and (130c) for λ i <λ j are automatically fulfilled by (cid:101) A ( z ) , which can be shown φ (1) φ (1) ξ η G ( ξ , η ) cb φ (1) ξ η G ( ξ , η ) φ (1) φ (1) ξ η G ( ξ , η ) cb φ (1) ξ η G ( ξ , η ) Fig. 3. Solution of the canonical kernel equations for µ c = 2 with thecolor proportional to the surface hight ranging from blue to yellow. The solidlines represent the BCs, whereas the dashed line at G ( ξ , η ) is theartificial BC with g f ( η ) = 0 needed for well-posedness. The red lines inthe spatial domains for the elements G ( ξ , η ) and G ( ξ , η ) arelines of separation. Furthermore, b = φ (1) − φ (1) and c = φ (1) + φ (1) specify the corresponding spatial domains. using the reciprocity relation between K ( z, ζ ) and L ( z, ζ ) (cf.[12, Ch. 4.5]). The remaining BVP has the same form as(24), so that there exists a piecewise C -solution L ( z, ζ ) byTheorem 2. Hence, (128) is bounded so that the proof ofthe theorem follows from Theorem 1 by utilizing standardarguments. VI. E XAMPLE
Consider a system of two coupled PIDEs (1) with the systemparameters Λ( z ) = (cid:20) − sin(2 πz ) 00 + z cos(2 πz ) (cid:21) , Φ( z ) = 0 A ( z ) = (cid:20) z + z (cid:21) , A ( z ) = (cid:20) z − zz − z (cid:21) F ( z, ζ ) = (cid:20) e z + ζ e z − ζ − e − ( z − ζ ) e − ( z + ζ ) (cid:21) (132a)and the differential operators θ [ x ( t )](0) = (cid:20) (cid:21) ∂ z x (0 , t ) + (cid:20) − (cid:21) x (0 , t ) (132b) θ [ x ( t )](1) = (cid:20) (cid:21) ∂ z x (1 , t ) + (cid:20) (cid:21) x (1 , t ) (132c)specifying the BCs, which is open-loop unstable. Note that λ > λ is valid for (132a). The left BC (132b) specifies oneDirichlet and one Robin BC on the unactuated boundary. Onthe other side, (132c) specifies Neumann actuation for the first EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 14 z ζ K ( z, ζ ) z ζ K ( z, ζ ) z ζ K ( z, ζ ) z ζ K ( z, ζ ) Fig. 4. Kernel elements in the original coordinates for µ c = 2 with thecolor proportional to the surface hight ranging from blue to yellow. The solidlines represent the BCs, whereas the dashed line at K ( z, ζ ) is the artificialBC needed for well-posedness. The red lines in the spatial domains for theelements K ( z, ζ ) and K ( z, ζ ) are lines of separation. state with a coupling of the second state which is actuatedby a Dirichlet BC. The controller is parametrized by µ c = 2 and the degree of freedom for specifying the artificial BC ischosen as g f ( η ) = 0 (see (76e)). For the target system (11a)decoupled BCs at both sides are imposed. In particular, thesame BC as in (132b) and ˜ θ [˜ x ( t )](1) = (cid:20) (cid:21) ∂ z ˜ x (1 , t ) + (cid:20) (cid:21) ˜ x (1 , t ) (133)hold.In order to solve the kernel equations, the successiveapproximation was numerically implemented in M ATLAB .Thereby, the spatial variables z and ζ were discretised using51 grid points each. The resulting grids for the different ( ξ, η ) -coordinate systems were resampled to ensure a minimum pointdistance for numerical performance. The successive approxima-tion is stopped as soon as max ξ,η,i,j max( | ∆ G lij | , | ∆ H lij | ) < − , which occurs after 11 iteration steps.Figure 3 shows the solutions G ij ( ξ ij , η ij ) , i, j = 1 , , ofthe canonical kernel equations (76) in the respective spatial do-mains. Those are triangular domains for the diagonal elements.For the off-diagonal elements, the lower boundary evolves inthe fourth quadrant because of the mutually different diffusioncoefficients and is a strictly monotonically decreasing non-straight line as the diffusion coefficients are spatially-varying.The BCs required for the mapping into the target system aremarked by the solid lines, whereas an artificial BC has to beassigned for the element G (dotted line). Since the spatialdomains of the elements G and G also cover the fourthquadrant, these kernel elements are only defined piecewise. . . − tz x ( z , t ) µ c = 2 . . − tz x ( z , t ) µ c = 2 . . ( µ max − t t (cid:107) x ( t ) (cid:107) () µ c = 2 . . ( µ max − t t (cid:107) x ( t ) (cid:107) () µ c = 8 Fig. 5. Solution of the closed-loop system for x ( z ) = (cid:2) sin( π z ) sin( πz ) − π cos( π z ) (cid:3) (cid:62) . The first row depicts the closed-loopstate profile for µ c = 2 , whereas the plots in the second row contain the L -norms of the solution for µ c = 2 and µ c = 8 with µ max = − . . This is indicated by the line of separation (red line), which isgiven by η = 0 and η = 0 . On these lines, the obtainedkernel elements are still continuous, but not differentiable.In Figure 4, the solution of the kernel equations K ij ( z, ζ ) , i, j = 1 , , in the original coordinates is depicted. Here, theline of separation (red line) is no longer a straight line, buta strictly monotonically increasing curve due to the nonlinearchange of coordinates.For the simulation, the plant was discretised using a finite-element method with 102 grid points for each state. Tovisualize the influence of the controller parameter µ c , Figure 5shows the closed-loop state profiles for µ c = 2 and the related L -norms for µ c = 2 and µ c = 8 . In both cases, the initialvalue is x ( z ) = (cid:2) sin( π z ) sin( πz ) − π cos( π z ) (cid:3) (cid:62) , whichis compatible to the BCs. The result verifies that exponentialstability is achieved and that a desired closed-loop stabilitymargin can be assigned by µ c .VII. C ONCLUDING REMARKS
An interesting research topic is to take a convective couplingin the PIDEs and a coupling at the unactuated boundaryinto account. As far as the solution of the kernel equationsis concerned such an extension seems to be directly possi-ble. However, since the aforementioned couplings have tobe included in the target system, this yields, in general, abidirectionally coupled system of parabolic PDEs. Hence,a different proof of well-posedness and stability is needed.For this, the corresponding results in [26] are of interest.In order to determine output feedback controllers the designof backstepping observers for the considered system class iscurrently under investigation.
EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 15 A PPENDIX
Proof of Theorem 1 . Introduce the operators A i h = λ i d z h − µ c h , µ c > µ max , i = 1 , , . . . , n , with D ( A i ) = { h ∈ L (0 , | h, d z h abs. cont. , d z h ∈ L (0 , , e (cid:62) i θ [ P (cid:62) e i h ](0) = e (cid:62) i ˜ θ [ P (cid:62) e i h ](1) = 0 } in theHilbert space H = L (0 , with the usual inner productand e i ∈ R n denoting the i -th unit vector. Therein thelargest eigenvalue µ max ∈ R of A i + µ c I , i = 1 , , . . . , n ,exists because − ( A i + µ c I ) are Sturm-Liouville operators (see[5]). With this, the operator A h = [ A h . . . A n h n ] (cid:62) , D ( A ) = D ( A ) ⊕ . . . ⊕ D ( A n ) can be defined. Furthermore,let ∆ h = − (cid:101) A ∗ ( R d z h (0) + R h (0)) with D (∆) = D ( A ) .Based on these preparations the target system (16) can beformulated in form of the abstract initial value problem (IVP) ˙˜ x ( t ) = ( A + ∆) ˜ x ( t ) , t > , ˜ x (0) = ˜ x ∈ D ( A ) , on thestate space X = ( L (0 , n with the L -inner product (cid:104)· , ·(cid:105) .Therein, A is the infinitesimal generator of an analytic C -semigroup , because the operators A i are the generators ofanalytic C -semigroups. The latter property follows from thefact that A i are Riesz-spectral operators satisfying the sectorcondition of [4, Ex. 2.18]. Moreover, the operator ∆ is A -bounded (see [11, Ch. IV, §
1, Sec. 2]) and thus also A -compact due to its finite-dimensional range (see [11, Rem. IV/1.13]).This implies an A -bound equal to zero (see [6, Lem. III/2.16]).Thus, the perturbed operator A + ∆ is the generator of ananalytic C -semigroup by [6, Th. III/2.10]. Hence, the IVP(16) is well-posed. The operators A i have a compact resolvent,since they have no eigenvalue at 0 (see [18, Th. 7.5.4]). Thus,also A has a compact resolvent, so that the A -boundednessof ∆ implies that A + ∆ has a compact resolvent (see [11,Th. IV/3.17]). Therefore, the spectrum of the operator A + ∆ is discrete (see [11, Th. III/6.29]). The spectrum determinedgrowth assumption (SDGA) holds for A + ∆ because it is agenerator of an analytic semigroup (see [23]). This and thetriangular structure of A + ∆ yield the decay rate µ c − µ max for the semigroup related to A + ∆ and thus the stability resultof the theorem. (cid:4) R EFERENCES[1] A. Baccoli, Y. Orlov, and A. Pisano, “On the boundary control of coupledreaction-diffusion equations having the same diffusivity parameters,”
Proc. CDC, Los Angeles, California, USA , pp. 5222–5228, 2014.[2] A. Baccoli and A. Pisano, “Anticollocated backstepping observer designfor a class of coupled reaction-diffusion PDEs,”
J. Control Science andEng. , vol. 2015, pp. 1–10, 2015.[3] A. Baccoli, A. Pisano, and Y. Orlov, “Boundary control of coupledreaction-diffusion processes with constant parameters,”
Automatica ,vol. 54, pp. 80–90, 2015.[4] R. F. Curtain and H. J. Zwart,
An Introduction to Infinite-DimensionalLinear Systems Theory . New York: Springer-Verlag, 1995.[5] C. Delattre, D. Dochain, and J. Winkin, “Sturm-Liouville systems areRiesz-spectral systems,”
Int. J. Appl. Math. Comput. Sci. , vol. 13, pp.481–484, 2003.[6] K.-J. Engel and R. Nagel,
One-parameter Semigroups for LinearEvolution Equations . New York: Springer-Verlag, 2000.[7] L. Hu, F. Di Meglio, R. Vazquez, and M. Krstic, “Control of homodi-rectional and general heterodirectional linear coupled hyperbolic PDEs,”
IEEE Trans. Autom. Control , vol. 61, pp. 3301–3314, 2016.[8] L. Hu, R. Vazquez, F. Di Meglio, and M. Krstic, “Boundary exponentialstabilization of 1-D inhomogeneous quasilinear hyperbolic systems,”
Submitted for publication, available at arXiv.org , 2015. [9] L. Jadachowski, T. Meurer, and A. Kugi, “Backstepping observers forlinear PDEs on higher-dimensional spatial domains,”
Automatica , vol. 51,pp. 85–97, 2015.[10] K. Jensen and W. Ray, “The bifurcation behavior of tubular reactors,”
Chem. Eng. Sci. , vol. 37, pp. 199–222, 1982.[11] T. Kato,
Perturbation Theory for Linear Operators . Berlin: Springer-Verlag, 1995.[12] M. Krstic and A. Smyshlyaev,
Boundary Control of PDEs — A Courseon Backstepping Designs . Philadelphia: SIAM, 2008.[13] B. Liu, D. Boutat, and D. Liu, “Backstepping observer-based outputfeedback control for a class of coupled parabolic PDEs with differentdiffusions,”
Syst. Control Lett. , vol. 97, pp. 61–69, 2016.[14] W. Liu, “Boundary feedback stabilization of an unstable heat equation,”
SIAM J. Control Optim. , vol. 42, pp. 1033–1042, 2003.[15] T. Meurer,
Control of Higher-Dimensional PDEs . Berlin: Springer-Verlag, 2013.[16] T. Meurer and A. Kugi, “Tracking control for boundary controlledparabolic PDEs with varying parameters: Combining backstepping anddifferential flatness,”
Automatica , vol. 45, pp. 1182–1194, 2009.[17] T. Myint-U and L. Debnath,
Linear Partial Differential Equations forScientists and Engineers . Boston: Birkh¨auser, 2007.[18] A. Naylor and G. Sell,
Linear Operator Theory in Engineering andScience . New York: Springer-Verlag, 1982.[19] Y. Orlov, A. Pisano, A. Pilloni, and E. Usai, “Output feedback stabiliza-tion of coupled reaction-diffusion processes with constant parameters,”
SIAM J. Control Optim. , vol. 55, pp. 4112–4155, 2017.[20] W. Ray,
Advanced Process Control . New York: McGraw-Hill, 1981.[21] A. Smyshlyaev and M. Krstic, “Closed-form boundary state feedbacksfor a class of 1-D partial integro-differential equations,”
IEEE Trans.Autom. Control , vol. 49, pp. 2185–2202, 2004.[22] ——, “On control design for PDEs with space-dependent diffusivity ortime-dependent reactivity,”
Automatica , vol. 41, pp. 1601–1608, 2005.[23] R. Triggiani, “On the stabilizability problem in Banach space,”
J. Math.Anal. Appl. , vol. 52, pp. 383–403, 1975.[24] R. Vazquez and M. Krstic, “Control of 1-D parabolic PDEs with Volterranonlinearities, Part I: Design,”
Automatica , vol. 44, pp. 2778–2790, 2008.[25] ——, “Control of 1-D parabolic PDEs with Volterra nonlinearities, PartII: Analysis,”
Automatica , vol. 44, pp. 2791–2803, 2008.[26] ——, “Boundary control of coupled reaction-advection-diffusion systemswith spatially-varying coefficients,”
IEEE Trans. Autom. Control , vol. 62,pp. 2026–2033, 2017.
Joachim Deutscher received the Dipl.-Ing. (FH)degree in electrical engineering from the Univer-sity of Applied Sciences W¨urzburg-Schweinfurt-Aschaffenburg, Germany, in 1996, the Dipl.-Ing.Univ. degree in electrical engineering, the Dr.-Ing.and the Dr.-Ing. habil. degrees both in automaticcontrol from the Friedrich-Alexander UniversityErlangen-Nuremberg (FAU), Germany, in 1999, 2003and 2010, respectively.From 2003–2010 he was a Senior Researcher atthe Chair of Automatic Control (FAU), in 2011 hewas appointed Associate Professor and since 2017 he is a Professor at thesame university. Currently, he is Head of the Infinite-Dimensional SystemsGroup at the Chair of Automatic Control (FAU).His research interests include control of distributed-parameter systems andcontrol theory for nonlinear lumped-parameter systems with applications inmechatronic systems.Dr. Deutscher has co-authored a book on state feedback control forlinear lumped-parameter systems: Design of Observer-Based Compensators(Springer, 2009) and is author of the book: State Feedback Control ofDistributed-Parameter Systems (in German) (Springer, 2012).
EEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. Y, JULY 2017 16