Feedback Stabilization of Non-Uniform Spatial ... - Semantic Scholar
2013 American Control Conference (ACC) Washington, DC, USA, June 17-19, 2013
Feedback Stabilization of Non-uniform Spatial Pattern in Reaction-Diffusion Systems Kenji Kashima∗ , Toshiyuki Ogawa† , Tatsunari Sakurai‡ Abstract— In this paper, we formulate and solve feedback stabilization problem of unstable non-uniform spatial pattern in reaction-diffusion systems. By considering spatial spectrum dynamics, we obtain a finite dimensional approximation that takes over the semi-passivity of the original partial differential equation. By virtue of this property, we can show the diffusive coupling in the spatial frequency domain achieves the desired pattern formation.
I. I NTRODUCTION The goal of this paper is to formulate and solve feedback stabilization problem of unstable non-uniform spatial pattern in reaction-diffusion systems [1], [2], [3], [4], [5]. See also [6], [7], [8], [9], [10], [11] for existing results related to control theory. We begin with a brief introduction of pattern generation mechanism in reaction-diffusion systems. Let spatiotemporal variable z(t, ξ) where the second variable ξ ∈ Ω is a spatial variable. Then, we consider the following partial differential equation ∂ z(t, ξ) = f (z(t, ξ)) + D∆z(t, ξ), ξ ∈ Ω (1) ∂t (with suitable boundary condition) with non-negative diagonal matrix D. We show an illustrative example in which z = [u, v]′ and Ω is a two-dimensional rectangular domain equipped with a periodic boundary condition. The origin of the reaction term (represented by f ) is globally asymptotically stable; see the following sections for more detail. Figure 1 is snapshots of u(t, ξ) of (1), where an initial state close to z eq ≡ 0 was randomly generated. The spatio-temporal pattern gradually goes away from the trivial equilibrium pattern, and finally converges to a roll pattern. On the other hand, we can see transient, but long lasting enough to observe, other specific patterns. In view of this, we attempt stabilize such patterns that cannot exist stably by spatially distributed feedback control. In this paper, this problem is formulated suitably, and then solved based on diffusive coupling in the spatial frequency domain. We briefly compare the contribution of this paper with existing results in controls community. Distributed parameter ∗ Graduate School of Engineering Science, Osaka University, Osaka, Japan [email protected] † Graduate School of Advanced Study of Mathematical Sciences, Meiji University, Kanagawa, Japan [email protected] ‡ Department of Physics, Graduate School of Science, Chiba University, Chiba, Japan [email protected] This research is partially supported by the Aihara Innovative Mathematical Modelling Project, the Japan Society for the Promotion of Science (JSPS) through the “Funding Program for World-Leading Innovative R&D on Science and Technology (FIRST Program),” initiated by the Council for Science and Technology Policy (CSTP).
978-1-4799-0176-0/$31.00 ©2013 AACC
system theory has investigated control of semilinear systems, including reaction-diffusion systems, from various aspects e.g., [12], [13], [14], [15]. However, the control objective is usually the stabilization of the trivial equilibrium z eq . This shows a clear contrast to our case where the exact profile of target pattern is unavailable. Though Arcak investigated non-uniform spatial pattern [6], [7], detailed discussion on the resulting pattern and control system design are not fully explored. The organization of this paper is as follows. Section II gives prerequisite for analysis of reaction-diffusion systems. In Section III, we formulate feedback control problem of spatial pattern formation. In Section IV, we reformulate and solve the problem based on spatial spectrum dynamics. Numerical examples are given in Section V to show the effectiveness of the obtained result. Section VI makes some concluding remarks. II. P RELIMINARY A. Reaction-diffusion systems In this paper, we consider the following form of reactiondiffusion system: the spatial domain is Ω := [0, Lx ] × [0, Ly ].
(2)
and spatio-temporal state variable is two-dimensional [ ] u(t, x, y) z(t, x, y) = ∈ R2 . v(t, x, y) Then, the spatio-temporal dynamics is given by ∂u = a11 u − a12 v − u3 + du ∆u ∂t ∂v = a21 u − a22 v + dv ∆v ∂t
(3) (4)
for all (x, y) ∈ Ω with the periodic boundary condition { z(t, 0, y) = z(t, Lx , y), for all (t, x, y). (5) z(t, x, 0) = z(t, x, Ly ), This reaction-diffusion system has a trivial equilibrium pattern z eq (x, y) ≡ 0 on Ω.
(6)
Even if the reaction term is globally stable, this trivial equilibrium of reaction-diffusion system (3), (4) is not necessarily stable.
3765
(a) t = 0
(f) t = 100
(b) t = 5
(g) t = 500
(c) t = 10
(h) t = 1000
(d) t = 20
(i) t = 1500
B. Spatial spectral analysis of uniform pattern instability
zm (t) :=
um (t) vm (t)
∫
A. Motivating example We take
∗
z(t, x, y)pm (x, y) dxdy ∈ C
:= Ω
pm (x, y) := e and
2
2πj( mLxxx +
z(t, x, y) =
2πj( mLxxx +
my y Ly )
(13) (14) (15)
.
m1 = (4, 0), m2 = (2, 2), m3 = (2, −2),
∗ zm = z−m for all m ∈ Z2
(8)
since the dynamics of our interest is real-valued. Then, it is useful to investigate the dynamics of each spatial component {zm (t)}m∈Z2 instead of {z(t, x, y)}(x,y)∈Ω . Let us consider the localized dynamics around the trivial pattern z eq . It should be emphasized that each wave number has its decoupled local dynamics: (9)
] , 0 dv
(10) ] ,
(16)
For these parameters, Am has one (real) unstable eigenvalue for m = ±mi
Note that
a ¯11 −a12 Am := A − sm D =: a21 −¯ a22 [ ] [ a11 −a12 du A := , D := a21 −a22 0 )2 ( )2 ( 2πmy 2πmx + . sm := Lx Ly
,
Lx 8π Lx = √ , Ly = √ . 1.8 3
m∈Z2
d zm (t) = Am zm (t), dt [
]
and
my y Ly )
zm (t)e
1 −1 3 −2
du = 0.2, dv = 1.5.
(7)
is its complex conjugate. This satisfies ∑
[ A=
for wave number m = (mx , my ) ∈ Z2 where
p∗m
(k) t = 3000
III. P ROBLEM F ORMULATION
In this section, the diffusion-induced local instability of the trivial equilibrium pattern is analyzed based on spatial spectral analysis. Consider the spatial Fourier transform ]
(j) t = 2000
Spatial pattern formation in reaction-diffusion systems: normalized snapshots of u(t, x, y).
Fig. 1.
[
(e) t = 50
(11) (12)
It should be emphasized that when du ̸= dv , stability of A does not necessarily guarantee stability of Am . In such a case, the corresponding spatial wave pm grows around the z eq . Further discussion on the pattern formation needs to consider the effect of nonlinearity.
(17)
and Am stable otherwise. The simulation in Fig. 1 was executed under these parameter settings. The initial patterns are randomly generated but sufficiently close to z eq . We can expect this reactiondiffusion system can generate three roll patterns corresponding to mi ’s in (17). Actually, all of observed patterns (including transient ones) look like superpositions of these roll patterns. B. Stabilization of unknown target pattern Having the example in the previous section in our mind, let us formulate stabilization problem of unstable spatial patterns. We define the set of wave numbers for which the local dynamics is unstable. Assumption 1: The finite set M ⊂ Z2 satisfies 1) Am has at least one eigenvalue in C+ if m ∈ M, 2) Am is stable if m ∈ / ±M := {±m : m ∈ M}, and 3) if m ∈ M, then −m ∈ / M. Because Am = A−m and (8), we imposed the condition 3) in order to avoid redundancy. It should be emphasized that sm1 = sm2 can hold for m1 ̸= m2 . Next, for feedback control problem, we assume that we can observe and also manipulate u in a spatially distributed
3766
manner. Thus, the controlled reaction-diffusion system is given by [ ] ∂ w z = f (z) + D∆z + , (18) 0 ∂t w(t, x, y) = W (u(t, x, y)) (19) Problem 1: Under definitions above and Assumption 1, find a feedback control law W (·) such that 1) z does not diverge, 2) zm for m ∈ / ±M asymptotically vanishes, 3) zm for m ∈ M converges to the same nonzero value, and 4) w(t, x, y) asymptotically vanishes. IV. F EEDBACK C ONTROL OF C ENTER M ANIFOLD DYNAMICS A. Reformulation on center manifold dynamics In view of Assumption 1-2), let us assume that zm is negligible for m ∈ / ±M to avoid infinite-dimensionality. We analyze the dynamics that zm should obey under the following technical assumption: Assumption 2: If n1 , n2 , n3 ∈ M satisfy n1 +n2 +n3 = m ∈ ±M, then at least one of ni ’s is equal to m. By putting zm ≈ 0 for m ∈ / ±M, ∑ u3 = un1 un2 un3 pn1 +n2 +n3 (x, y).
1’) zm is ultimately bounded, that is, there exists (initial state independent) R > 0 such that lim sup ∥zm ∥ < R for all m ∈ M, 3’-a) lim |zm (t) − zn (t)| = 0 for all m, n ∈ M
t→∞
(n1 , n2 , n3 ) = (m, n, −n) and its permutation, where n ∈ ±M is arbitrary. Therefore, we obtain the following approximation in the spatial frequency domain: ] [ ] [ d um um = Am (20) vm dt vm ∑ [ ] −um 3|um |2 + 6 |un |2 + wm + 0 n̸=m 0 ∫
where wm (t) :=
w(t, x, y)p∗m (x, y)dxdy.
Ω
It emphasized that, in the rest of this paper, zm = [ should ]be ′ um , vm represents the solution of finite dimensionally approximated system (20), instead of the spatial spectral component of the solution u, v to the original partial differential equation (3), (4). Now, Problem 1 is naturally rephrased in this domain: Problem 2: Consider the complex-valued system (20). Then, find a feedback control law
lim wm (t) = 0 for all m ∈ M. (23) For this problem, we can simply obtain a desired control law by considering the specific structure of Problem 2. Theorem 1: For Problem 2, consider the following diffusive coupling on M ∑ γmn (un (t) − um (t)) (24) wm (t) = σ t→∞
n∈M,n̸=m
where γmn ≥ 0. If γmn is associated to a strongly connected graph on M, then there exists σ > 0 such that the control law (24) with any strength σ > σ satisfies the following: • •
the conditions 1’) in Problem 2 hold. if sm = s¯ (see (12) for the definition) for all m ∈ M, then 3’-a), 3’-b), 4’) are also satisfied.
B. Proof of Theorem 1 The proof is similar to that of [16], [17]. However, we need to pay attention to the complex-valued dynamics and also the extra interaction term arising from u3 . For notational simplicity, denote zM = (zm )m∈M , uM = (um )m∈M , vM = (vm )m∈M u fm (uM , vm ) = a ¯11 um − a12 vm − um βm (uM ), v fm (uM , vm ) = a21 um − a ¯22 vm , ∑ |un |2 ≥ 0. βm (uM ) := 3|um |2 + 6 n̸=m, n∈M
1) Ultimate boundedness: In this section, we prove the first claim. First, we begin with semi-passivity like property for subsystems. Lemma 1: Define the real-valued functions defined on C2 ( ) a12 2 1 |u|2 + |v| (25) V (u, v) = 2 a21 and a12 a ¯22 2 |v| . a21
(26)
d ∗ V (um (t), vm (t)) ≤ Re[wm um ] − H(um , vm ) dt
(27)
H(u, v) = |u|2 (3|u|2 − a ¯11 ) + Then,
wm (t) = Wm (uM (t)) such that,
(22)
3’-b) the origin is locally unstable, and 4’)
n1 ,n2 ,n3 ∈±M
It follows from Assumption 2 and the orthogonality of {pm }m∈Z2 that the wave number m component appears only from the combination
(21)
t→∞
holds for all m ∈ M. 3767
Proof: By direct computation
2) Consensus: Next, we show the second claim. Define z˜i = zmi − zm1 , i = 2, · · · , |M|, z˜ := [˜ z2 , · · · , z˜|M| ]′
d V (um (t), vm (t)) dt ( [ u ]) a12 ∗ fm + wm ∗ = Re [um , v ] v fm a21 m ∗ = Re[wm um ] + |um |2 (¯ a11 − βm (uM )) − ∗ ≤ Re[wm um ] − H(um , vm ).
(and also u ˜i , v˜i , w ˜i and u ˜, v˜, w, ˜ similarly). It suffices to show z˜ converges to the origin. For • = u, v, define
a12 a ¯22 |vm |2 a21
• f˜• (uM , vM ) := [f˜2• (uM , vM ), · · · , f˜|M| (uM , vM )]′ (34) • • f˜i• (uM , vM ) := fm (uM , vmi ) − fm (uM , vm1 ). (35) i 1
Lemma 2: Take positive numbers νm such that ∑ νn γn,m = 0 for all m
(28)
n∈M
with
∑
γmm := −
γmn .
Lemma 3: For any matrix P ∈ R(|M|−1)×(|M|−1) such that ∥P ∥ ≤ 1, there exists a positive constant c3 such that Re(˜ u∗ P f˜u (uM , vM )) ≤ a12 ∥˜ u∥∥˜ v ∥ + c3 ∥˜ u∥2 (1 + ∥uM ∥2 ). (36)
(29)
Proof: We have
n̸=m
f˜iu =¯ a11 u ˜i − a12 v˜i − umi βmi (uM ) + um1 βm1 (uM ) ∑ =¯ a11 u ˜i − a12 v˜i − 6˜ ui |um |2
Then, ∑
νm Re[u∗m wm ] ≤ 0 for all m.
(30)
m∈M
+ 3(umi |umi |2 − um1 |um1 |2 ).
m∈M
Proof: The existence of νm follows from the PerronFrobenius theorem. We have ∑ ∗ νm (u∗m wm + wm um )
m
≤
∑ m
=
∑
νm (
m
∑
∑
u ˜i (|umi |2 + |um1 |2 ) + u ˜i u∗mi u∗m1 = umi |umi |2 − um1 |um1 |2 . Therefore, there exists c3 > 0 such that
γmn (|un |2 − |um |2 )
|f˜iu + a12 v˜i | ≤ c3 |˜ ui |(1 + ∥uM ∥2 ), i = 2, · · · , |M|
)
νn γnm
|um |2 = 0
and consequently
n∈M
∥f˜u ∥ ≤ a12 ∥˜ v ∥ + c3 ∥˜ u∥(1 + ∥uM ∥2 ).
where we used |um | + |un | − 2
For the last term,
n̸=m
2
(32) (33)
u∗m un
−
u∗n um
This readily gives the desired result. We are in the position to prove the consensus formation, that is, the convergence of z˜ to 0. By definition,
= |um − un | ≥ 0. 2
Hence, (30) holds. We are ready to show the boundedness. Define ∑ W (zM ) := νm V (um , vm ).
d u ˜ = f˜u (uM , vM ) + w. ˜ dt d v˜ = f˜v (uM , vM ). dt We can always take Hurwitz matrix L such that
m∈M
Then, by Lemma 1 and 2, ∑ d W (zM ) ≤ − νm H(zm ). dt
w ˜ = σL˜ u. (31)
m∈M
This inequality combined with the positivity of νm indicate d that we can take c such that dt W (zM ) is strictly negative outside of B2 := {zM : W (zM ) ≤ c}.
Take µ > 0 and a positive definite matrix P such that ∥P ∥ ≤ 1 satisfing the Lyapunov equation LP + P L′ = −2µI, so that
This means that B2 is positively invariant, and all solutions enter this compact set in finite time. This completes the proof of the ultimate boundedness. 3768
Re(˜ u∗ P L˜ u) = −µ∥˜ u∥2 .
Let us define V3 (˜ u, v˜) :=
) 1( ∗ u ˜ Pu ˜ + ∥¯ v ∥2 . 2
0.025
By direct computation, we have v˜i∗ (f v (umi , vmi ) − f v (um1 , vm1 )) = −¯ a22 |˜ vi |2 + a21 v˜i∗ u ˜i .
0.02
Therefore, by the ultimate boundedness, there exists c4 > 0 such that d V3 (˜ u(t), v˜(t)) dt ∑ ∗ u 2 ∗ = Re u ˜ P (f˜ (˜ u, v˜) + σL˜ u) + (−¯ a22 |˜ vi | + a21 v˜i u ˜ i )
0.015
0.01
0.005
i̸=1
( ∗ u ) = Re u ˜ P f¯ (˜ u, v˜) + a21 v˜∗ u ˜ −a ¯22 ∥˜ v ∥2 − σµ∥˜ u∥2
0 0
≤a12 ∥˜ u∥∥˜ v ∥ + c4 ∥˜ u∥2 − σµ∥˜ u∥2 − a ¯22 ∥˜ v ∥2 + a21 ∥˜ u∥˜ v∥
1000
2000
3000
4000
5000
6000
≤(c4 − σµ)∥˜ u∥2 + (a12 + a21 )∥˜ u∥∥˜ v∥ − a ¯22 ∥¯ v ∥2 . Fig. 2.
The right hand side is negative definite if c4 − σµ < 0, −4¯ a22 (c4 − σµ) > (a12 + a21 )2 . This means if 1 σ> µ
Time response of um : Uncontrolled.
0.025
) ( (a12 + a21 )2 c4 + =: σ, 4¯ a22
0.02
0.015
then V3 (˜ u, v˜) converge to 0 as t → ∞. Thus, we have 3’-a). Let A¯ := As¯. Then, the linearized dynamics of zM ¯ around the origin is given by
0.01
0.005
d zM = (I ⊗ A¯ + Γ ⊗ D1 )zM , dt [ ] 1 0 D1 = 0 0 { γmi ,mj i ̸= j (Γ)i,j = γmi ,mi in (29), i = j.
(37) 0 0
20
40
60
80
100
(38) (39)
Note Γ has 0 as an eigenvalue. This readily means ¯ ⊂ eig(I ⊗ A¯ + Γ ⊗ D1 ). eig(A) Recall that A¯ is unstable by Assumption 1-2). This completes the proof of 3’-b). Trivially, we have 4’) since w(t) also converges to 0 by definition. V. N UMERICAL EXAMPLE We consider system (20) arising from the motivating example in Section III. The initial state is generated from the spatial spectrum of the initial pattern in Fig. 1. Fig. 2 shows the time response without any control input wm ≡ 0, where • real- and imaginary-parts are shown by solid and dashed lines, and • unstable wave numbers m = (4, 0), (2, 2), (2, −2) are shown by blue, red, and magenta, respectively. As expected from the convergence to the roll pattern in Fig. 1, two spatial spectrum asymptotically vanishes. On the other hand, Figs. 3 and 4 show the time response for the controlled case with σ = 0.1. We can see all um (m ∈ M) converges to the same nonzero value, which corresponds
Fig. 3. Time response of um : Controlled by the proposed diffusive coupling in the spatial spectrum domain (t = 0 ∼ 100).
to the hexagonal pattern in the original reaction-diffusion system. Note that we cannot generate square pattern by partial diffusive coupling of the corresponding spatial spectrum, for example, by taking −1 1 0 Γ = 1 −1 0 . 0 0 0 In this case, the coupled spatial spectrum achieve consensus, however, possibly converges to 0, as shown in Fig. 5. Then, we observe a roll pattern as in Fig. 1. VI. C ONCLUSION In this paper, we formulated a feedback stabilization problem of unstable non-uniform spatial pattern in the reaction-diffusion systems. This problem was solved in the finite dimensionally approximated system. The proposed law, which is a diffusive coupling in the spatial spectrum, suitably achieves spectrum consensus and instability preservation. As a next step, we need to embed the results for the proposed control law onto the original partial differential
3769
0.025
0.02
0.015
0.01
0.005
0 0
500
1000
1500
Fig. 4. Time response of um : Controlled by the proposed diffusive coupling in the spatial spectrum domain (t = 0 ∼ 1500). 0.025
0.02
0.015
0.01
0.005
0 0
500
1000
1500
2000
2500
3000
Fig. 5. Time response of um : Controlled by the partial diffusive coupling in the spatial spectrum domain.
equation. The corresponding input pattern should be given as w(t, x, y) := ∑ 2σ Re pm m∈M
∑
(40)
[3] R. Kapral, Chemical Waves and Patterns, K. Showalter, Ed. Kluwer Academic Publishers, 1995. [4] J. D. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications. Springer, 2003. [5] A. S. Mikhailov and K. Showalter, “Control of waves, patterns and turbulence in chemical systems,” Physics Reports, vol. 425, no. 2-3, pp. 79–194, Mar. 2006. [6] M. Arcak, “Certifying spatially uniform behavior in reaction-diffusion PDE and compartmental ODE systems,” Automatica, vol. 47, no. 6, pp. 1219–1229, Jun. 2011. [7] J. Hsia, W. A. Holtz, D. C. Huang, M. Arcak, and M. M. Maharbiz, “A quenched oscillator network for pattern formation in gene expression,” in American Control Conference (ACC), 2011, 2011, pp. 2284–2289. [8] S. Luther, F. H. Fenton, B. G. Kornreich, A. Squires, P. Bittihn, D. Hornung, M. Zabel, J. Flanders, A. Gladuli, L. Campoy, E. M. Cherry, G. Luther, G. Hasenfuss, V. I. Krinsky, A. Pumir, R. F. Gilmour, and E. Bodenschatz, “Low-energy control of electrical turbulence in the heart.” Nature, vol. 475, no. 7355, pp. 235–9, Jul. 2011. [9] K. Konishi, M. Takeuchi, and T. Shimizu, “Design of external forces for eliminating traveling wave in a piecewise linear FitzHugh-Nagumo model.” Chaos (Woodbury, N.Y.), vol. 21, no. 2, p. 023101, Jun. 2011. [10] L. M. Mu˜noz, J. F. Stockton, and N. F. Otani, “Applications of control theory to the dynamics and propagation of cardiac action potentials.” Annals of biomedical engineering, vol. 38, no. 9, pp. 2865–76, Sep. 2010. [11] C. Vilas, M. R. Garc´ıa, J. R. Banga, and A. A. Alonso, “Robust feed-back control of travelling waves in a class of reaction-diffusion distributed biological systems,” Physica D: Nonlinear Phenomena, vol. 237, no. 18, pp. 2353–2364, Sep. 2008. [12] Y. Smagina and M. Sheintuch, “Using Lyapunov’s direct method for wave suppression in reactive systems,” Systems & Control Letters, vol. 55, no. 7, pp. 566–572, Jul. 2006. [13] E. Fridman and Y. Orlov, “An LMI approach to boundary control of semilinear parabolic and hyperbolic systems,” Automatica, vol. 45, no. 9, pp. 2060–2066, Sep. 2009. [14] R. F. Curtain and H. J. Zwart, An Introduction to Infinite-Dimensional Linear Systems Theory. Springer-Verlag, 1995. [15] M. Krstic and A. Smyshlyaev, Boundary Control of PDEs: A Course on Backstepping Designs. Society for Industrial and Applied Mathematics, 2008. [16] E. Steur, I. Tyukin, and H. Nijmeijer, “Semi-passivity and synchronization of diffusively coupled neuronal oscillators,” Physica D: Nonlinear Phenomena, vol. 238, no. 21, pp. 2119–2128, Nov. 2009. [17] E. Steur and H. Nijmeijer, “Synchronization in networks of diffusively time-delay coupled (semi-)passive systems,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 58, no. 6, pp. 1358–1371, Jun. 2011.
γmn (un − um )
n̸=m, n∈M
where um for m ∈ M is defined by (7). We have already verified numerically that this control law achieves the expected pattern formation. Ongoing works contain • theoretical evaluation of the control effect in the original reaction-diffusion systems, • generalization of generated pattern, e.g., square pattern for the example above, and • generalization of the dynamics such as chemotaxis equations, Rayleigh-B´enard convection. R EFERENCES [1] A. M. Turing, “The Chemical Basis of Morphogenesis,” Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 237, no. 641, pp. 37–72, Aug. 1952. [2] A. T. Winfree, When Time Breaks Down. Princeton Univ. Press, 1987.
3770
Feedback Stabilization of Non-uniform Spatial Pattern in Reaction-Diffusion Systems Kenji Kashima∗ , Toshiyuki Ogawa† , Tatsunari Sakurai‡ Abstract— In this paper, we formulate and solve feedback stabilization problem of unstable non-uniform spatial pattern in reaction-diffusion systems. By considering spatial spectrum dynamics, we obtain a finite dimensional approximation that takes over the semi-passivity of the original partial differential equation. By virtue of this property, we can show the diffusive coupling in the spatial frequency domain achieves the desired pattern formation.
I. I NTRODUCTION The goal of this paper is to formulate and solve feedback stabilization problem of unstable non-uniform spatial pattern in reaction-diffusion systems [1], [2], [3], [4], [5]. See also [6], [7], [8], [9], [10], [11] for existing results related to control theory. We begin with a brief introduction of pattern generation mechanism in reaction-diffusion systems. Let spatiotemporal variable z(t, ξ) where the second variable ξ ∈ Ω is a spatial variable. Then, we consider the following partial differential equation ∂ z(t, ξ) = f (z(t, ξ)) + D∆z(t, ξ), ξ ∈ Ω (1) ∂t (with suitable boundary condition) with non-negative diagonal matrix D. We show an illustrative example in which z = [u, v]′ and Ω is a two-dimensional rectangular domain equipped with a periodic boundary condition. The origin of the reaction term (represented by f ) is globally asymptotically stable; see the following sections for more detail. Figure 1 is snapshots of u(t, ξ) of (1), where an initial state close to z eq ≡ 0 was randomly generated. The spatio-temporal pattern gradually goes away from the trivial equilibrium pattern, and finally converges to a roll pattern. On the other hand, we can see transient, but long lasting enough to observe, other specific patterns. In view of this, we attempt stabilize such patterns that cannot exist stably by spatially distributed feedback control. In this paper, this problem is formulated suitably, and then solved based on diffusive coupling in the spatial frequency domain. We briefly compare the contribution of this paper with existing results in controls community. Distributed parameter ∗ Graduate School of Engineering Science, Osaka University, Osaka, Japan [email protected] † Graduate School of Advanced Study of Mathematical Sciences, Meiji University, Kanagawa, Japan [email protected] ‡ Department of Physics, Graduate School of Science, Chiba University, Chiba, Japan [email protected] This research is partially supported by the Aihara Innovative Mathematical Modelling Project, the Japan Society for the Promotion of Science (JSPS) through the “Funding Program for World-Leading Innovative R&D on Science and Technology (FIRST Program),” initiated by the Council for Science and Technology Policy (CSTP).
978-1-4799-0176-0/$31.00 ©2013 AACC
system theory has investigated control of semilinear systems, including reaction-diffusion systems, from various aspects e.g., [12], [13], [14], [15]. However, the control objective is usually the stabilization of the trivial equilibrium z eq . This shows a clear contrast to our case where the exact profile of target pattern is unavailable. Though Arcak investigated non-uniform spatial pattern [6], [7], detailed discussion on the resulting pattern and control system design are not fully explored. The organization of this paper is as follows. Section II gives prerequisite for analysis of reaction-diffusion systems. In Section III, we formulate feedback control problem of spatial pattern formation. In Section IV, we reformulate and solve the problem based on spatial spectrum dynamics. Numerical examples are given in Section V to show the effectiveness of the obtained result. Section VI makes some concluding remarks. II. P RELIMINARY A. Reaction-diffusion systems In this paper, we consider the following form of reactiondiffusion system: the spatial domain is Ω := [0, Lx ] × [0, Ly ].
(2)
and spatio-temporal state variable is two-dimensional [ ] u(t, x, y) z(t, x, y) = ∈ R2 . v(t, x, y) Then, the spatio-temporal dynamics is given by ∂u = a11 u − a12 v − u3 + du ∆u ∂t ∂v = a21 u − a22 v + dv ∆v ∂t
(3) (4)
for all (x, y) ∈ Ω with the periodic boundary condition { z(t, 0, y) = z(t, Lx , y), for all (t, x, y). (5) z(t, x, 0) = z(t, x, Ly ), This reaction-diffusion system has a trivial equilibrium pattern z eq (x, y) ≡ 0 on Ω.
(6)
Even if the reaction term is globally stable, this trivial equilibrium of reaction-diffusion system (3), (4) is not necessarily stable.
3765
(a) t = 0
(f) t = 100
(b) t = 5
(g) t = 500
(c) t = 10
(h) t = 1000
(d) t = 20
(i) t = 1500
B. Spatial spectral analysis of uniform pattern instability
zm (t) :=
um (t) vm (t)
∫
A. Motivating example We take
∗
z(t, x, y)pm (x, y) dxdy ∈ C
:= Ω
pm (x, y) := e and
2
2πj( mLxxx +
z(t, x, y) =
2πj( mLxxx +
my y Ly )
(13) (14) (15)
.
m1 = (4, 0), m2 = (2, 2), m3 = (2, −2),
∗ zm = z−m for all m ∈ Z2
(8)
since the dynamics of our interest is real-valued. Then, it is useful to investigate the dynamics of each spatial component {zm (t)}m∈Z2 instead of {z(t, x, y)}(x,y)∈Ω . Let us consider the localized dynamics around the trivial pattern z eq . It should be emphasized that each wave number has its decoupled local dynamics: (9)
] , 0 dv
(10) ] ,
(16)
For these parameters, Am has one (real) unstable eigenvalue for m = ±mi
Note that
a ¯11 −a12 Am := A − sm D =: a21 −¯ a22 [ ] [ a11 −a12 du A := , D := a21 −a22 0 )2 ( )2 ( 2πmy 2πmx + . sm := Lx Ly
,
Lx 8π Lx = √ , Ly = √ . 1.8 3
m∈Z2
d zm (t) = Am zm (t), dt [
]
and
my y Ly )
zm (t)e
1 −1 3 −2
du = 0.2, dv = 1.5.
(7)
is its complex conjugate. This satisfies ∑
[ A=
for wave number m = (mx , my ) ∈ Z2 where
p∗m
(k) t = 3000
III. P ROBLEM F ORMULATION
In this section, the diffusion-induced local instability of the trivial equilibrium pattern is analyzed based on spatial spectral analysis. Consider the spatial Fourier transform ]
(j) t = 2000
Spatial pattern formation in reaction-diffusion systems: normalized snapshots of u(t, x, y).
Fig. 1.
[
(e) t = 50
(11) (12)
It should be emphasized that when du ̸= dv , stability of A does not necessarily guarantee stability of Am . In such a case, the corresponding spatial wave pm grows around the z eq . Further discussion on the pattern formation needs to consider the effect of nonlinearity.
(17)
and Am stable otherwise. The simulation in Fig. 1 was executed under these parameter settings. The initial patterns are randomly generated but sufficiently close to z eq . We can expect this reactiondiffusion system can generate three roll patterns corresponding to mi ’s in (17). Actually, all of observed patterns (including transient ones) look like superpositions of these roll patterns. B. Stabilization of unknown target pattern Having the example in the previous section in our mind, let us formulate stabilization problem of unstable spatial patterns. We define the set of wave numbers for which the local dynamics is unstable. Assumption 1: The finite set M ⊂ Z2 satisfies 1) Am has at least one eigenvalue in C+ if m ∈ M, 2) Am is stable if m ∈ / ±M := {±m : m ∈ M}, and 3) if m ∈ M, then −m ∈ / M. Because Am = A−m and (8), we imposed the condition 3) in order to avoid redundancy. It should be emphasized that sm1 = sm2 can hold for m1 ̸= m2 . Next, for feedback control problem, we assume that we can observe and also manipulate u in a spatially distributed
3766
manner. Thus, the controlled reaction-diffusion system is given by [ ] ∂ w z = f (z) + D∆z + , (18) 0 ∂t w(t, x, y) = W (u(t, x, y)) (19) Problem 1: Under definitions above and Assumption 1, find a feedback control law W (·) such that 1) z does not diverge, 2) zm for m ∈ / ±M asymptotically vanishes, 3) zm for m ∈ M converges to the same nonzero value, and 4) w(t, x, y) asymptotically vanishes. IV. F EEDBACK C ONTROL OF C ENTER M ANIFOLD DYNAMICS A. Reformulation on center manifold dynamics In view of Assumption 1-2), let us assume that zm is negligible for m ∈ / ±M to avoid infinite-dimensionality. We analyze the dynamics that zm should obey under the following technical assumption: Assumption 2: If n1 , n2 , n3 ∈ M satisfy n1 +n2 +n3 = m ∈ ±M, then at least one of ni ’s is equal to m. By putting zm ≈ 0 for m ∈ / ±M, ∑ u3 = un1 un2 un3 pn1 +n2 +n3 (x, y).
1’) zm is ultimately bounded, that is, there exists (initial state independent) R > 0 such that lim sup ∥zm ∥ < R for all m ∈ M, 3’-a) lim |zm (t) − zn (t)| = 0 for all m, n ∈ M
t→∞
(n1 , n2 , n3 ) = (m, n, −n) and its permutation, where n ∈ ±M is arbitrary. Therefore, we obtain the following approximation in the spatial frequency domain: ] [ ] [ d um um = Am (20) vm dt vm ∑ [ ] −um 3|um |2 + 6 |un |2 + wm + 0 n̸=m 0 ∫
where wm (t) :=
w(t, x, y)p∗m (x, y)dxdy.
Ω
It emphasized that, in the rest of this paper, zm = [ should ]be ′ um , vm represents the solution of finite dimensionally approximated system (20), instead of the spatial spectral component of the solution u, v to the original partial differential equation (3), (4). Now, Problem 1 is naturally rephrased in this domain: Problem 2: Consider the complex-valued system (20). Then, find a feedback control law
lim wm (t) = 0 for all m ∈ M. (23) For this problem, we can simply obtain a desired control law by considering the specific structure of Problem 2. Theorem 1: For Problem 2, consider the following diffusive coupling on M ∑ γmn (un (t) − um (t)) (24) wm (t) = σ t→∞
n∈M,n̸=m
where γmn ≥ 0. If γmn is associated to a strongly connected graph on M, then there exists σ > 0 such that the control law (24) with any strength σ > σ satisfies the following: • •
the conditions 1’) in Problem 2 hold. if sm = s¯ (see (12) for the definition) for all m ∈ M, then 3’-a), 3’-b), 4’) are also satisfied.
B. Proof of Theorem 1 The proof is similar to that of [16], [17]. However, we need to pay attention to the complex-valued dynamics and also the extra interaction term arising from u3 . For notational simplicity, denote zM = (zm )m∈M , uM = (um )m∈M , vM = (vm )m∈M u fm (uM , vm ) = a ¯11 um − a12 vm − um βm (uM ), v fm (uM , vm ) = a21 um − a ¯22 vm , ∑ |un |2 ≥ 0. βm (uM ) := 3|um |2 + 6 n̸=m, n∈M
1) Ultimate boundedness: In this section, we prove the first claim. First, we begin with semi-passivity like property for subsystems. Lemma 1: Define the real-valued functions defined on C2 ( ) a12 2 1 |u|2 + |v| (25) V (u, v) = 2 a21 and a12 a ¯22 2 |v| . a21
(26)
d ∗ V (um (t), vm (t)) ≤ Re[wm um ] − H(um , vm ) dt
(27)
H(u, v) = |u|2 (3|u|2 − a ¯11 ) + Then,
wm (t) = Wm (uM (t)) such that,
(22)
3’-b) the origin is locally unstable, and 4’)
n1 ,n2 ,n3 ∈±M
It follows from Assumption 2 and the orthogonality of {pm }m∈Z2 that the wave number m component appears only from the combination
(21)
t→∞
holds for all m ∈ M. 3767
Proof: By direct computation
2) Consensus: Next, we show the second claim. Define z˜i = zmi − zm1 , i = 2, · · · , |M|, z˜ := [˜ z2 , · · · , z˜|M| ]′
d V (um (t), vm (t)) dt ( [ u ]) a12 ∗ fm + wm ∗ = Re [um , v ] v fm a21 m ∗ = Re[wm um ] + |um |2 (¯ a11 − βm (uM )) − ∗ ≤ Re[wm um ] − H(um , vm ).
(and also u ˜i , v˜i , w ˜i and u ˜, v˜, w, ˜ similarly). It suffices to show z˜ converges to the origin. For • = u, v, define
a12 a ¯22 |vm |2 a21
• f˜• (uM , vM ) := [f˜2• (uM , vM ), · · · , f˜|M| (uM , vM )]′ (34) • • f˜i• (uM , vM ) := fm (uM , vmi ) − fm (uM , vm1 ). (35) i 1
Lemma 2: Take positive numbers νm such that ∑ νn γn,m = 0 for all m
(28)
n∈M
with
∑
γmm := −
γmn .
Lemma 3: For any matrix P ∈ R(|M|−1)×(|M|−1) such that ∥P ∥ ≤ 1, there exists a positive constant c3 such that Re(˜ u∗ P f˜u (uM , vM )) ≤ a12 ∥˜ u∥∥˜ v ∥ + c3 ∥˜ u∥2 (1 + ∥uM ∥2 ). (36)
(29)
Proof: We have
n̸=m
f˜iu =¯ a11 u ˜i − a12 v˜i − umi βmi (uM ) + um1 βm1 (uM ) ∑ =¯ a11 u ˜i − a12 v˜i − 6˜ ui |um |2
Then, ∑
νm Re[u∗m wm ] ≤ 0 for all m.
(30)
m∈M
+ 3(umi |umi |2 − um1 |um1 |2 ).
m∈M
Proof: The existence of νm follows from the PerronFrobenius theorem. We have ∑ ∗ νm (u∗m wm + wm um )
m
≤
∑ m
=
∑
νm (
m
∑
∑
u ˜i (|umi |2 + |um1 |2 ) + u ˜i u∗mi u∗m1 = umi |umi |2 − um1 |um1 |2 . Therefore, there exists c3 > 0 such that
γmn (|un |2 − |um |2 )
|f˜iu + a12 v˜i | ≤ c3 |˜ ui |(1 + ∥uM ∥2 ), i = 2, · · · , |M|
)
νn γnm
|um |2 = 0
and consequently
n∈M
∥f˜u ∥ ≤ a12 ∥˜ v ∥ + c3 ∥˜ u∥(1 + ∥uM ∥2 ).
where we used |um | + |un | − 2
For the last term,
n̸=m
2
(32) (33)
u∗m un
−
u∗n um
This readily gives the desired result. We are in the position to prove the consensus formation, that is, the convergence of z˜ to 0. By definition,
= |um − un | ≥ 0. 2
Hence, (30) holds. We are ready to show the boundedness. Define ∑ W (zM ) := νm V (um , vm ).
d u ˜ = f˜u (uM , vM ) + w. ˜ dt d v˜ = f˜v (uM , vM ). dt We can always take Hurwitz matrix L such that
m∈M
Then, by Lemma 1 and 2, ∑ d W (zM ) ≤ − νm H(zm ). dt
w ˜ = σL˜ u. (31)
m∈M
This inequality combined with the positivity of νm indicate d that we can take c such that dt W (zM ) is strictly negative outside of B2 := {zM : W (zM ) ≤ c}.
Take µ > 0 and a positive definite matrix P such that ∥P ∥ ≤ 1 satisfing the Lyapunov equation LP + P L′ = −2µI, so that
This means that B2 is positively invariant, and all solutions enter this compact set in finite time. This completes the proof of the ultimate boundedness. 3768
Re(˜ u∗ P L˜ u) = −µ∥˜ u∥2 .
Let us define V3 (˜ u, v˜) :=
) 1( ∗ u ˜ Pu ˜ + ∥¯ v ∥2 . 2
0.025
By direct computation, we have v˜i∗ (f v (umi , vmi ) − f v (um1 , vm1 )) = −¯ a22 |˜ vi |2 + a21 v˜i∗ u ˜i .
0.02
Therefore, by the ultimate boundedness, there exists c4 > 0 such that d V3 (˜ u(t), v˜(t)) dt ∑ ∗ u 2 ∗ = Re u ˜ P (f˜ (˜ u, v˜) + σL˜ u) + (−¯ a22 |˜ vi | + a21 v˜i u ˜ i )
0.015
0.01
0.005
i̸=1
( ∗ u ) = Re u ˜ P f¯ (˜ u, v˜) + a21 v˜∗ u ˜ −a ¯22 ∥˜ v ∥2 − σµ∥˜ u∥2
0 0
≤a12 ∥˜ u∥∥˜ v ∥ + c4 ∥˜ u∥2 − σµ∥˜ u∥2 − a ¯22 ∥˜ v ∥2 + a21 ∥˜ u∥˜ v∥
1000
2000
3000
4000
5000
6000
≤(c4 − σµ)∥˜ u∥2 + (a12 + a21 )∥˜ u∥∥˜ v∥ − a ¯22 ∥¯ v ∥2 . Fig. 2.
The right hand side is negative definite if c4 − σµ < 0, −4¯ a22 (c4 − σµ) > (a12 + a21 )2 . This means if 1 σ> µ
Time response of um : Uncontrolled.
0.025
) ( (a12 + a21 )2 c4 + =: σ, 4¯ a22
0.02
0.015
then V3 (˜ u, v˜) converge to 0 as t → ∞. Thus, we have 3’-a). Let A¯ := As¯. Then, the linearized dynamics of zM ¯ around the origin is given by
0.01
0.005
d zM = (I ⊗ A¯ + Γ ⊗ D1 )zM , dt [ ] 1 0 D1 = 0 0 { γmi ,mj i ̸= j (Γ)i,j = γmi ,mi in (29), i = j.
(37) 0 0
20
40
60
80
100
(38) (39)
Note Γ has 0 as an eigenvalue. This readily means ¯ ⊂ eig(I ⊗ A¯ + Γ ⊗ D1 ). eig(A) Recall that A¯ is unstable by Assumption 1-2). This completes the proof of 3’-b). Trivially, we have 4’) since w(t) also converges to 0 by definition. V. N UMERICAL EXAMPLE We consider system (20) arising from the motivating example in Section III. The initial state is generated from the spatial spectrum of the initial pattern in Fig. 1. Fig. 2 shows the time response without any control input wm ≡ 0, where • real- and imaginary-parts are shown by solid and dashed lines, and • unstable wave numbers m = (4, 0), (2, 2), (2, −2) are shown by blue, red, and magenta, respectively. As expected from the convergence to the roll pattern in Fig. 1, two spatial spectrum asymptotically vanishes. On the other hand, Figs. 3 and 4 show the time response for the controlled case with σ = 0.1. We can see all um (m ∈ M) converges to the same nonzero value, which corresponds
Fig. 3. Time response of um : Controlled by the proposed diffusive coupling in the spatial spectrum domain (t = 0 ∼ 100).
to the hexagonal pattern in the original reaction-diffusion system. Note that we cannot generate square pattern by partial diffusive coupling of the corresponding spatial spectrum, for example, by taking −1 1 0 Γ = 1 −1 0 . 0 0 0 In this case, the coupled spatial spectrum achieve consensus, however, possibly converges to 0, as shown in Fig. 5. Then, we observe a roll pattern as in Fig. 1. VI. C ONCLUSION In this paper, we formulated a feedback stabilization problem of unstable non-uniform spatial pattern in the reaction-diffusion systems. This problem was solved in the finite dimensionally approximated system. The proposed law, which is a diffusive coupling in the spatial spectrum, suitably achieves spectrum consensus and instability preservation. As a next step, we need to embed the results for the proposed control law onto the original partial differential
3769
0.025
0.02
0.015
0.01
0.005
0 0
500
1000
1500
Fig. 4. Time response of um : Controlled by the proposed diffusive coupling in the spatial spectrum domain (t = 0 ∼ 1500). 0.025
0.02
0.015
0.01
0.005
0 0
500
1000
1500
2000
2500
3000
Fig. 5. Time response of um : Controlled by the partial diffusive coupling in the spatial spectrum domain.
equation. The corresponding input pattern should be given as w(t, x, y) := ∑ 2σ Re pm m∈M
∑
(40)
[3] R. Kapral, Chemical Waves and Patterns, K. Showalter, Ed. Kluwer Academic Publishers, 1995. [4] J. D. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications. Springer, 2003. [5] A. S. Mikhailov and K. Showalter, “Control of waves, patterns and turbulence in chemical systems,” Physics Reports, vol. 425, no. 2-3, pp. 79–194, Mar. 2006. [6] M. Arcak, “Certifying spatially uniform behavior in reaction-diffusion PDE and compartmental ODE systems,” Automatica, vol. 47, no. 6, pp. 1219–1229, Jun. 2011. [7] J. Hsia, W. A. Holtz, D. C. Huang, M. Arcak, and M. M. Maharbiz, “A quenched oscillator network for pattern formation in gene expression,” in American Control Conference (ACC), 2011, 2011, pp. 2284–2289. [8] S. Luther, F. H. Fenton, B. G. Kornreich, A. Squires, P. Bittihn, D. Hornung, M. Zabel, J. Flanders, A. Gladuli, L. Campoy, E. M. Cherry, G. Luther, G. Hasenfuss, V. I. Krinsky, A. Pumir, R. F. Gilmour, and E. Bodenschatz, “Low-energy control of electrical turbulence in the heart.” Nature, vol. 475, no. 7355, pp. 235–9, Jul. 2011. [9] K. Konishi, M. Takeuchi, and T. Shimizu, “Design of external forces for eliminating traveling wave in a piecewise linear FitzHugh-Nagumo model.” Chaos (Woodbury, N.Y.), vol. 21, no. 2, p. 023101, Jun. 2011. [10] L. M. Mu˜noz, J. F. Stockton, and N. F. Otani, “Applications of control theory to the dynamics and propagation of cardiac action potentials.” Annals of biomedical engineering, vol. 38, no. 9, pp. 2865–76, Sep. 2010. [11] C. Vilas, M. R. Garc´ıa, J. R. Banga, and A. A. Alonso, “Robust feed-back control of travelling waves in a class of reaction-diffusion distributed biological systems,” Physica D: Nonlinear Phenomena, vol. 237, no. 18, pp. 2353–2364, Sep. 2008. [12] Y. Smagina and M. Sheintuch, “Using Lyapunov’s direct method for wave suppression in reactive systems,” Systems & Control Letters, vol. 55, no. 7, pp. 566–572, Jul. 2006. [13] E. Fridman and Y. Orlov, “An LMI approach to boundary control of semilinear parabolic and hyperbolic systems,” Automatica, vol. 45, no. 9, pp. 2060–2066, Sep. 2009. [14] R. F. Curtain and H. J. Zwart, An Introduction to Infinite-Dimensional Linear Systems Theory. Springer-Verlag, 1995. [15] M. Krstic and A. Smyshlyaev, Boundary Control of PDEs: A Course on Backstepping Designs. Society for Industrial and Applied Mathematics, 2008. [16] E. Steur, I. Tyukin, and H. Nijmeijer, “Semi-passivity and synchronization of diffusively coupled neuronal oscillators,” Physica D: Nonlinear Phenomena, vol. 238, no. 21, pp. 2119–2128, Nov. 2009. [17] E. Steur and H. Nijmeijer, “Synchronization in networks of diffusively time-delay coupled (semi-)passive systems,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 58, no. 6, pp. 1358–1371, Jun. 2011.
γmn (un − um )
n̸=m, n∈M
where um for m ∈ M is defined by (7). We have already verified numerically that this control law achieves the expected pattern formation. Ongoing works contain • theoretical evaluation of the control effect in the original reaction-diffusion systems, • generalization of generated pattern, e.g., square pattern for the example above, and • generalization of the dynamics such as chemotaxis equations, Rayleigh-B´enard convection. R EFERENCES [1] A. M. Turing, “The Chemical Basis of Morphogenesis,” Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 237, no. 641, pp. 37–72, Aug. 1952. [2] A. T. Winfree, When Time Breaks Down. Princeton Univ. Press, 1987.
3770
