Feedback Stabilization for a coupled PDE-ODE Production System
Vanessa Baumgärtner Simone Göttlich Stephan Knapp








[1]University of Mannheim, Department of Mathematics, 68131 Mannheim, Germany (goettlich@uni-mannheim.de).



Feedback Stabilization for a coupled PDE-ODE Production System
    Vanessa Baumgärtner[1],   Simone  Göttlich[1],   Stephan Knapp[1]
    March 26, 2019
=====================================================================




We consider an interlinked production model consisting of conservation laws (PDE) coupled to ordinary differential equations (ODE).
Our focus is the analysis of control laws for the coupled system and corresponding stabilization questions of equilibrium dynamics in the 
presence of disturbances. These investigations are carried out using an appropriate Lyapunov function on the theoretical and numerical level. 
The discrete L^2-stabilization technique allows to derive a mixed feedback law that is able to ensure exponential stability 
also in bottleneck situations. All results are accompanied by computational examples.


AMS Classification: 65Mxx, 93D05, 90B30
 



Keywords: Feedback stabilization, Lyapunov function, coupled PDE-ODE system




§ INTRODUCTION

Mathematical models to describe production systems have gained a lot of attention
during the past few decades. In particular, fluid-like models based on scalar hyperbolic conservation laws (PDE)
have been developed and analyzed, see for example 
<cit.>. 
The extension to networks has been originally introduced in <cit.>, wherein, as coupling conditions, queues in terms of
ordinary differential equations (ODE) have been installed. This approach immediately leads to a coupled PDE-ODE dynamics, see <cit.> 
for more details. 

In this article, we are concerned with the stabilization of such a coupled PDE-ODE system.
The dynamics within each processing unit is given by the linear advection equation with positive velocity while
queues in front of the unit are introduced to avoid congestion. Therefore, the queues measure the difference between the in-and outflow
into a processing unit and therefore obey a nonlinear ordinary differential equation.  
In recent years, the stability of linear coupled systems <cit.> has been investigated intensively.
However, there is only little literature available on the stability analysis for nonlinear coupled systems, see for example <cit.>.
Therefore, our purpose is to focus on the stability analysis of the coupled PDE-ODE production system from a theoretical as well as
discrete (or numerical) point of view. 
Due to the fact that the coupled PDE-ODE system can be interpreted as a coupled boundary value problem,
stabilization results for hyperbolic equations can be applied and extended, see <cit.>.
Therein, analytical results for sufficiently smooth solutions in connection with boundary control have been obtained. 
The underlying tools for the study of those problems are Lyapunov functions stabilizing the deviation from steady states in suitable norms, e.g. L^2 in space. Exponential decay of a continuous Lyapunov function under a so-called dissipative boundary condition has been proven in <cit.>.
Also, explicit decay rates for numerical schemes have only recently been established. In <cit.>, exponential decay on a finite time horizon has been developed for non-conservative (first-order) schemes and in <cit.> for discretizations of linear systems only.  
Different to the already existing results for linear feedback laws, we derive a mixed feedback law based on the numerical discretization, 
that is capable to tackle instances, where queues are still filled and hence the ODE affects the stability behavior.


This paper is organized as follows: in Section <ref>, we introduce the production model that presented in <cit.> and provide a suitable numerical discretization. We make use of the theory of Lyapunov functions (Section <ref>) to derive a feedback law that ensures the exponential stability for the coupled PDE-ODE model in Section <ref>. The theoretical results are accompanied by various simulation results. 
Since we are also interested in the exponential stability of the numerical solution, we prove its exponential stability and additionally provide decay rates in Section <ref>. On the discrete level, we obtain a so-called mixed feedback law to deal with queueing situations, i.e. non-empty queue loads.




§ MODEL EQUATIONS AND NUMERICAL DISCRETIZATION


First, we recall the production model taken from <cit.>. For simplicity, we stick to the special case of a serial system without any branches. The model considers a queue in front of every processing unit e∈{1,...,m}. Furthermore, the maximal capacity μ_e>0 and the processing velocity v_e>0 of each unit are constant parameters. 
The dynamics of the product density ρ_e(t,x) defined on a segment [a_e,b_e] follows <cit.>:

	
    ∂_tρ_e(t,x)+∂_xf(ρ_e(t,x))=0,

	
    f_e(ρ_e(t,x))=min{v_eρ_e(t,x),μ_e},    ∀ x∈[a_e,b_e],   t∈ℝ_+,

	
    ρ_e(0,x)=ρ_e,0(x),


where f_e denotes the flux function. We define the inflow g_in,e at x=a_e in front of processing unit e as the outflow of the predecessor:

    g_in,e(t) f_e-1(ρ_e-1(t,b_e-1)), for e=2,...,m.

In the case e=1, an externally given inflow profile g_in,1(t)=f_in(t) is assumed.
Due to the possibility of different maximal capacities μ_e, it might occur that the inflow 
can not be processed immediately. Therefore, we introduce a time-dependent function q_e(t) describing the load of the queue at processing unit e.
Each queue satisfies the following ordinary differential equation:

    ∂_tq_e(t)=g_in,e(t)-g_out,e(t), with g_out,e(t)  f_e(ρ_e(t,a_e))

describing the outflow from queue e to processor e at time t, cf. Figure <ref>. 


Obviously, the outgoing flux g_out,e depends on the queue load. 
So if the queue is empty, the outflow is the minimum of the incoming flux g_in,e and the maximum processing capacity μ_e. In the first case, the queue remains empty, whereas in the second case, the queue starts to increase. In the non-empty case, the outgoing flux is equal to the maximal capacity. Summarizing, we get

    g_out,e(t) f_e(ρ_e(t,a_e))=min{f_e-1(ρ_e-1(t,b_e-1)),μ_e},   q_e(t)=0,
    μ_e,   q_e(t)>0.

Equations (<ref>) and (<ref>) ensure that the production model is well-defined, provided that f_in(t) and ρ_e,0(x) are of total bounded variation, see <cit.>.

To avoid discontinuities in the derivative of the queue length, a smoothed out version of (<ref>) has been derived in <cit.> and reads as

    f_e(ρ_e(t,a_e))≈min{μ_e,q_e(t)/ϵ}, with ϵ≪1.

Provided v_eρ_e,0≤μ_e, (<ref>) reduces to a linear advection equation

    ∂_tρ_e+v_e∂_xρ_e=0.
 
For the numerical investigations later on, we introduce a first-order discretization of the coupled production model. 
From now on, we assume that x∈[0,l], where l is the uniform length of all processors. 
It is advantageous to rewrite the model in terms of fluxes only to apply the Lyapunov stability concepts in Section <ref>.
By defining 

    Λ=diag(v_1,…,v_m),    f(t,x)=(f_1(t,x),…,f_m(t,x))^T ,   ∂_tq(t)=(∂_tq_1(t),…,∂_tq_m(t))^T,
     
     g_in(t)=(g_in,1(t),…,g_in,m(t))^T  and    g_out(t)=(g_out,1(t),…,g_out,m(t))^T,

we are able to write the model in the following compact form:

	
    ∂_tf(t,x)+Λ∂_xf(t,x)=0, x∈[0,l],  t∈ℝ_+,

	
    ∂_tq(t)=g_in(t)-g_out(t),

	

equipped with initial values f_e,0(x)=v_eρ_e,0(x), q(0)=q_0≥0. Note that (<ref>) is already given in characteristic form, see <cit.>.
 
For the discretization of (<ref>) we apply the left-sided upwind method due to the strictly positive velocities v_e.
The equation for the queues is approximated by an explicit Euler scheme as proposed in <cit.>. 
We use an equidistant grid with constant time step size τ and constant space step size h. To describe the boundary values, we use one ghost cell at the left. We consider N+1 cells with cell centers x_j=(j+1/2)h for j∈{-1,…,N-1}. The interfaces of the cells are located at x_j-1/2=jh for j∈{-1,...,N-1}. The considered time domain is [0,T], whereby the discrete time is denoted by t_k=kτ for k∈{0,...,K-1} such that (K-1)τ=T. Moreover, the CFL condition is assumed to be fulfilled, i.e.,

    v̅τ/h≤ 1,   v̅max_e∈{1,...,m}v_e.


For further investigations, we assume that 
q_1^k=0 ∀ k, which is satisfied if f_in^k≤μ_1 ∀ k. 
Then, the discretized equations for j∈{0,...,N-1} and k∈{0,...,K-2} read:

	
    f_e,j^k+1   =f_e,j^k-τ/hv_e(f_e,j^k-f_e,j-1^k),   for e∈{1,...,m},
    
    	
    	
    	
    	f_e,-1^k   =g_out,e^k=min{f_e-1,N-1^k,μ_e},   q_e^k=0,
    μ_e,                        q_e^k>0,
    	   for e≠1,
    
    	f_1,-1^k   =g_out,1^k=f_in^k,
    
    	q_e^k+1   =q_e^k+τ(g_in,e^k-g_out,e^k), with   for e∈{1,...,m},
    g_in,e^k   =f_e-1,N-1^k, e≠1, 
    
    	f_in^k,         else,
    
    	f_e,j^0    = f_0,e, q_e^0= q_0,e,   for e∈{1,...,m}.

	
 
Here, (<ref>),(<ref>) are the discretized boundary conditions that are determined by the queue or the control input. 
We note that we distinguish between the source node (e=1) and internal nodes (e = 2,…,m). Due to assumption f_in^k≤μ_1 ∀ k, 
it is valid to set g_out,1(t) equal to the inflow profile f_in(t). The discretized production network model is referred to as PNM in the following.





§ FEEDBACK STABILIZATION FOR THE CONTINUOUS SYSTEM



We now investigate the stabilization of the continuous model (<ref>). In this section, starting with the Lyapunov stability, we discuss results for a serial network consisting of only two processors and provide numerical illustrations. We remark that the extension to more processing units is
straightforward.


 §.§ Lyapunov Stability


In the following, we analyze the exponential stability of the continuous system (<ref>). 
In contrast to <cit.>, we are faced with queues and therefore investigate the exponential stability in the norm

    ‖(f(t,·),q(t))‖=√(‖ f(t,·)‖_L^2((0,l);ℝ^m)^2+|q(t)|^2),

where |·| is the Euclidean norm. The introduced norm (<ref>) is the natural choice of the product space L^2((0,l);ℝ^m)×ℝ^m motivated by <cit.>.

	The system (<ref>)-(<ref>) is exponentially stable in the sense of the norm (<ref>), if there exist ν>0 and C>0 such that for every initial condition f_0∈ L^2((0,l);ℝ^m) and q(0)∈ℝ^m the solution to (<ref>) satisfies
	
    ‖ f(t,·)‖_L^2((0,l);ℝ^m)^2+|q(t)|^2≤ Ce^-ν t(‖ f(0,·)‖_L^2((0,l);ℝ^m)^2+|q(0)|^2), ∀ t∈ℝ_+.
	
	


Note that we are particularly interested in the stabilization of the trivial steady state f≡0 and q≡0. Nevertheless, the model also allows for non-trivial states, which read as
f ≡ z, q ≡ 0, where

    z=(z_1, …,z_m), with z_e=ζmin_e(μ_e) for 0≤ζ<1.

By translation f̃=f-z, we recover results for the non-trivial steady states from the trivial steady state again.
For the feedback analysis we make use of a Lyapunov function approach. Since we consider a serial system with queues, the Lyapunov function already known for 
hyperbolic equations (see e.g. <cit.>) needs to be adapted. We define the following Lyapunov function candidate:

    V(t)   =V_1(t)+V_2(t)
       =∫_0^l f^T(t,x)P(x)f(t,x) dx+q^T(t)Q(t)q(t),
 with the weighting matrices
P(x)=diag(p_1e^-η _1x,...,p_me^-η_mx), p_e>0,η_e>0 and
Q(t)=diag(c_1e^-η̃_1 v_1t,…,c_me^-η̃_mv_mt), c_e>0,η̃_e>0. In order to include the ODEs for the queues, we add the quadratic form V_2(t)=q(t)^TQ(t)q(t). If not stated otherwise, we assume η_e=η and η̃_e=η̃, independent of the individual properties. 


For the discrete stability analysis in Section <ref>, we define the discrete version of the Lyapunov function in (<ref>) as

    V^k=∑_j=0^N-1∑_e=1^m(f_e,j^k)^2p_ee^-η_e x_e,jh+∑_e=1^m(q_e^k)^2c_ee^-η̃_ev_et_k,   with η_e,η̃_e,p_e,c_e>0.






 §.§ Analytical Feedback Stabilization


The goal is to derive a feedback law using the Lyapunov function from (<ref>). 
Figure <ref> illustrates the key idea of such a closed-loop system, where
a percentage of the outflow f_2(t,l) is fed back into the system.
We aim to derive a control input u_1(t) that depends on the outflow of the system such that exponential stability is ensured. 
We call this control input a feedback law.




By considering a control input u_1(t) as an inflow into the system, equations (<ref>) and (<ref>) of the numerical discretization are adapted as follows:


    f_1,-1^k   =g_out,1^k=u_1^k,
    
    q_e^k+1   =q_e^k+τ(g_in,e^k-g_out,e^k), with   for e∈{1,...,m},
    g_in,e^k   =f_e-1,N-1^k, e≠1,
     
    u_1^k,          else.




Since so far, only the non-queue case has been studied in the literature <cit.>. However, 
difficulties arise when queues start to grow up. Therefore, we need to adapt the already established results accordingly.
The following lemma presents a condition that ensures the exponential decay of the Lyapunov function (<ref>).

	The continuous solution satisfies
	
    V̇≤-ν V,

	i.e., exponential decay of V is obtained, if 
	
    -f(t,l)^TΛ P(l)f(t,l)+f(t,0)^TΛ P(0)f(t,0)+2∂_tq(t)^TQ(t)q(t)≤0

	is fulfilled ∀ t ∈ℝ_+. 
	The decay rate ν is given by
	
    ν=min(η,η̃)min_e(v_e)>0.

	


	In the following, we assume that f and q are sufficiently smooth. We obtain:
	
    [                                                                              V̇                                                                        (<ref>)=                       ∫_0^l(∂_tf^TP(x)f+f^TP(x)∂_tf)dx+2∂_tq^TQ(t)q+q^T∂_tQ(t)q;                                                                           [1ex]                                                                        (<ref>)=                  ∫_0^l(-∂_xf^TΛ P(x)f-f^TP(x)Λ∂_xf)dx+2∂_tq^TQ(t)q-η̃q^TΛ Q(t)q;                                                                           [1ex]                                                                               =            -∫_0^l∂_x[f^TΛ P(x)f]dx-η∫_0^lf^TΛ P(x)fdx+2∂_tq^TQ(t)q-η̃q^TΛ Q(t)q;                                                                           [1ex]                                                                               ≤ -[f^TΛ P(x)f]_0^l-ηmin_e(v_e)∫_0^lf^T P(x)fdx+2∂_tq^TQ(t)q-η̃min_e(v_e)q^TQ(t)q;                                                                           [1ex]                                                                               ≤                         -min(η,η̃)min_e(v_e)V(t)-[f^TΛ P(x)f]_0^l+2∂_tq^TQ(t)q. ]

	To obtain the exponential decay of V, we deduce the condition
	
    -[f^TΛ P(x)f]_0^l+2∂_tq^TQ(t)q≤0.





Next, we derive a feedback law by using Lemma <ref> and by imposing the feedback control u_1(t) at processor 1.
For the serial network in Figure <ref> we have to show

    -v_1f_1(t,l)^2p_1e^-η l+v_1f_1(t,0)^2p_1e^-η0-v_2f_2(t,l)^2p_2e^-η l+v_2f_2(t,0)^2p_2e^-η 0
       +2q_1(t)(g_in,1(t)-g_out,1(t))c_1e^-η̃v_1t+2q_2(t)(g_in,2(t)-g_out,2(t))c_2e^-η̃v_2t≤0.

Due to the assumption q_1(t)=0, the outflow into processor 1 is equal to the control input u_1(t), i.e., f_1(t,0)=u_1(t). This leads to

    v_1u_1(t)^2p_1   ≤v_1f_1(t,l)^2p_1e^-η l+v_2f_2(t,l)^2p_2e^-η l
        -v_2f_2(t,0)^2p_2-2q_2(t)g_in,2(t)c_2e^-η̃v_2t+2q_2(t)g_out,2(t)c_2e^-η̃v_2t.

In the next step, we set p_1=p_2=c_1=c_2=1 and use the approximate version of g_out,2(t)=f_2(t,0) in (<ref>) for ϵ≪1 to ensure enough regularity:

    u_1(t)^2   ≤ f_1(t,l)^2e^-η l+v_2/v_1f_2(t,l)^2e^-η l-v_2/v_1min{q_2(t)/ϵ,μ_2}^2
       -21/v_1q_2(t)f_1(t,l)e^-η̃v_2t+21/v_1q_2(t)min{q_2(t)/ϵ,μ_2}e^-η̃v_2t X(t).

Thus, the control input u_1(t) is bounded from above by X(t). This is only valid if X(t) is non-negative. In the steady state f≡0 and q≡0, the control u_1=0 is obtained.
Assuming v_e=v for e∈{1,2}, we get:

    X(t)   = f_1(t,l)^2e^-η l+f_2(t,l)^2e^-η l-min{q_2(t)/ϵ,μ_2}^2
       -21/vq_2(t)f_1(t,l)e^-η̃vt+21/vq_2(t)min{q_2(t)/ϵ,μ_2}e^-η̃vt.


Assuming v_e=v and p_e=p=c_e=c=1, this result can be extended to a serial system with m processors:  

    X(t)   =∑_e=1^mf_e(t,l)^2e^-η l-∑_e=2^mmin{q_e(t)/ϵ,μ_e}^2-∑_e=2^m21/vq_e(t)f_e-1(t,l)e^-η̃vt
       +∑_e=2^m21/vq_e(t)min{q_e(t)/ϵ,μ_e}e^-η̃vt.


If X(t)<0, we may set u_1(t)=0. Otherwise, we set u_1(t)=√(X(t)) motivated by (<ref>). We observe that u_1(t) depends on the values of f at the boundaries and on the queue load. 
This leads to a linear feedback control, i.e., a percentage of the outflow that is fed back into the system,

    u_1(t) = κ f_m(t,l),

for κ>0. 	

In the case of vanishing queues, i.e. q_e(t)=0 for all t∈ℝ_+, the results from Lemma <ref> and (<ref>)–(<ref>) coincide with <cit.> by a direct comparison.  




 §.§ Simulation Results

In the following, we provide simulation results 
for the linear feedback law for queueing situations. 
If not stated otherwise, we use a time horizon of T=50 and h=τ=0.01 for the step sizes.
All queues shall be empty at t=0, i.e. q_e(0)=0.
We set v_e=1 and the length l of each processor is 1.


We consider the linear feedback control law (<ref>) for three processors and two active queues. 
Our first approach to generate queues is to set the initial condition for the fluxes equal to the maximal capacity.
We choose μ_1=10, μ_2=9, μ_3=8 and the initial conditions f_1(0,x)=10, f_2(0,x)=9, f_3(0,x)=8. The latter choice leads to queues at the second and third processor. The parameter for the weighting matrices of the Lyapunov function are η_e=η̃_e=0.1 and p_e=c_e=1. 
We denote the upper bound from the proof of Lemma <ref> by:

    V_up(t) V(0)exp(-min(η,η̃)min_e(v_e)t).

Figure <ref> shows that the estimated Lyapunov function for κ=0.1 has a kink and is even slightly increasing for a short time interval because of the queues in front of processor 2 and 3. 
The Lyapunov function for the queue (V_2) is illustrated in Figure <ref>. 



Zooming in, we observe that this Lyapunov function decreases before it starts to increase again. This is due to the chosen initial conditions that cause queue 3 to continue growing while queue 2 already shrinks. Even though the Lyapunov function for the flux is decreasing or temporarily constant, the increase in queue 3 can not be compensated. Furthermore, we can see that the Lyapunov function for the flux is temporarily constant during two time intervals which is caused by the existence of two non-empty queues. 
Summarizing, there exist scenarios using the linear feedback law, where we are not able to ensure asymptotic stability, i.e. a decreasing Lyapunov function, for all time steps. However, such kinks can be avoided applying a so-called mixed feedback law, see Section <ref>.




§ FEEDBACK STABILIZATION FOR THE DISCRETIZED MODEL


In this section, we study how the results from our theoretical investigations can be transferred to the
discretized equations, similar to <cit.>.
We present decay rates for the exponential stability and derive a new type of feedback law
that is able to avoid increasing Lyapunov functions in case of non-empty queues.

Considering a discrete Lyapunov function (<ref>) also requires the definition of a discrete norm. 
Therefore, we use for the stability analysis of the disctrized equations

    ‖(f,q)‖_h=√(‖ f‖_h,2^2+|q|^2),     (f,q)∈ℝ^m× N×ℝ^m,

where ‖ f‖^2_h,2=∑_j=0^N-1∑_e=1^m(f_e,j)^2h is the discrete L^2-norm of f.

Since we are interested in a feedback stabilization for the numerical discretization, we equip (<ref>) with

	
    g_in(t)=Gf(t,l),


where the matrix G∈ℝ^m× m is the inflow matrix with constant entries, coupling the inflow with the outflow of the system. 
The following theorem states under which conditions the numerical discretization for the serial network model is exponentially stable.

	
	Let T>0 be arbitrarily large but fixed, f_in^k ≤μ_1 
	and the feedback law f_in^k=u_1^k=f_1,-1^k satisfies for k∈{0,...,K-1} the inequality:  
	
    ∑_e=1^mv_e((f_e,-1^k)^2p_eexp(-η_e x_e,0)-(f_e,N-1^k)^2p_eexp(-η_ex_e,N-1))
       +∑_e=1^m(2q_e^k(g_in,e^k-g_out,e^k)+τ(g_in,e^k-g_out,e^k)^2)c_eexp(-η̃_ev_et_k)exp(-η̃_ev_eτ)≤0.

	Then, provided the CFL condition (<ref>) is fulfilled, the following holds: there exist η_e>0 and η̃_e>0 such that the numerical scheme (<ref>)-(<ref>) is exponentially stable, i.e., the numerical solution satisfies
	
    V^k+1≤exp(-ντ(k+1))V^0, for k∈{0,…,K-2},

	where the decay rate ν is given by:
	
    ν=min{1/hmin_e(v_e(1-exp(-η_eh))),1/hmax_e(v_e)min_e(1-exp(-η̃_ev_eτ))}>0.
 
	Furthermore, f_e,j^k+1 and q_e^k+1 are exponentially stable in the sense of the discrete norm in (<ref>), i.e., 
	
    ∑_j=0^N-1∑_e=1^m(f_e,j^k+1)^2h+∑_e=1^m(q_e^k+1)^2≤ Cexp(-ντ(k+1))(∑_j=0^N-1∑_e=1^m(f_e,j^0)^2h+∑_e=1^m(q_e^0)^2),

	for k∈{0,...,K-2} and C>0.

 

	In order to obtain the exponential stability of the numerical solution, we show that the discrete time derivative of the discrete Lyapunov function in (<ref>) fulfills:
	
    V^k+1-V^k/τ≤-ν V^k, for a decay rate ν>0.

	By inserting the discrete Lyapunov function at time t^k and t^k+1, we obtain:
	
	
    V^k+1-V^k/τ   =∑_j=0^N-1∑_e=1^m((f_e,j^k+1)^2-(f_e,j^k)^2)p_eexp(-η_ex_e,j)h/τ
       +∑_e=1^m((q_e^k+1)^2c_eexp(-η̃_ev_et_k+1)-(q_e^k)^2c_eexp(-η̃_ev_et_k))1/τ
       =C_1+C_2,

	with
	
    C_1   =∑_j=0^N-1∑_e=1^m((f_e,j^k+1)^2-(f_e,j^k)^2)p_eexp(-η_ex_e,j)h/τ,
    
    	C_2   =∑_e=1^m((q_e^k+1)^2c_eexp(-η̃_ev_et_k+1)-(q_e^k)^2c_eexp(-η̃_ev_et_k))1/τ.

	At first, we have a closer look at C_1 by using the numerical discretization scheme (<ref>)-(<ref>) for the advection equation:
	
    [                                                                                     C_1                                                                                (<ref>)= ∑_j=0^N-1∑_e=1^m((f_e,j^k-τ/hv_e(f_e,j^k-f_e,j-1^k))^2-(f_e,j^k)^2)p_eexp(-η_ex_e,j)h/τ;                                                                                                                                                                               =               ∑_j=0^N-1∑_e=1^m((1-τ/hv_e)f_e,j^k+τ/hv_ef_e,j-1^k)^2p_eexp(-η_ex_e,j)h/τ;                                                                                                                                                                               -                                        ∑_j=0^N-1∑_e=1^m(f_e,j^k)^2p_eexp(-η_ex_e,j)h/τ. ]

	Applying the CFL condition, 
	we obtain
	
    C_1    ≤∑_j=0^N-1∑_e=1^m((1-τ/hv_e)(f_e,j^k)^2+τ/hv_e(f_e,j-1^k)^2)p_eexp(-η_ex_e,j)h/τ
        - ∑_j=0^N-1∑_e=1^m(f_e,j^k)^2p_eexp(-η_ex_e,j)h/τ
        = ∑_j=0^N-1∑_e=1^m(τ/hv_e(f_e,j-1^k)^2-τ/hv_e(f_e,j^k)^2)p_eexp(-η_ex_e,j)h/τ
        = ∑_e=1^mv_e∑_j=0^N-1((f_e,j-1^k)^2-(f_e,j^k)^2)p_eexp(-η_ex_e,j).

	In the next step, we make use of an index shift that requires the knowledge of the cell center x_N outside of the domain which is not needed in our discretization, since we only consider positive velocities. Therefore, by assuming x_e,N=x_e,N-1, we get:
	
    [                                                                                        C_1                                                                                          ≤ ∑_e=1^mv_e(∑_j=-1^N-2(f_e,j^k)^2p_eexp(-η_ex_e,j+1)-∑_j=0^N-1(f_e,j^k)^2p_eexp(-η_ex_e,j));                                                                                                                                                                                     =                    ∑_e=1^mv_e(∑_j=0^N-1(f_e,j^k)^2(p_eexp(-η_ex_e,j+1)-p_eexp(-η_ex_e,j)));                                                                                                                                                                                     +                   ∑_e=1^mv_e((f_e,-1^k)^2p_eexp(-η_ex_e,0)-(f_e,N-1^k)^2p_eexp(-η_ex_e,N));                                                                                                                                                                        x_e,N=x_e,N-1=                          ∑_e=1^mv_e(∑_j=0^N-1(exp(-η_eh)-1)(f_e,j^k)^2(p_eexp(-η_ex_e,j)));                                                                                                                                                                                     +                ∑_e=1^mv_e((f_e,-1^k)^2p_eexp(-η_ex_e,0)-(f_e,N-1^k)^2p_eexp(-η_ex_e,N-1)). ]

	Therefore, we have
	
    C_1≤ S_1+S_2,

	with
	
    S_1   =∑_e=1^mv_e(∑_j=0^N-1(exp(-η_eh)-1)(f_e,j^k)^2p_eexp(-η_ex_e,j)),
    
    	S_2   =∑_e=1^mv_e((f_e,-1^k)^2p_eexp(-η_e x_e,0)-(f_e,N-1^k)^2p_eexp(-η_ex_e,N-1)).

Regarding C_2, using the numerical discretization, or more specifically, the explicit Euler method for the queues, we obtain:
	
    [                                                                                            C_2                                                                                       (<ref>)= ∑_e=1^m((q_e^k+τ(g_in,e^k-g_out,e^k))^2c_eexp(-η̃_ev_et_k+1)-(q_e^k)^2c_eexp(-η̃_ev_et_k))1/τ.;                                                                                                ]

	In the next step, we use the binomial formula, which yields:
	
    [                                                                                                C_2                                                                                                  =                ∑_e=1^m((q_e^k)^2c_eexp(-η̃_ev_et_k)exp(-η̃_ev_eτ)-(q_e^k)^2c_eexp(-η̃_ev_et_k))1/τ;                                                                                                                                                                                                     + ∑_e=1^m(2q_e^kτ(g_in,e^k-g_out,e^k)+τ^2(g_in,e^k-g_out,e^k)^2)1/τc_eexp(-η̃_ev_et_k)exp(-η̃_ev_eτ);                                                                                                                                                                                                     =                                                                                           Z_1+Z_2, ]

	with
	
    Z_1   =∑_e=1^m(q_e^k)^2c_eexp(-η̃_ev_et_k)(exp(-η̃_ev_eτ)-1)1/τ,
    
    	Z_2   =∑_e=1^m(2q_e^k(g_in,e^k-g_out,e^k)+τ(g_in,e^k-g_out,e^k)^2)c_eexp(-η̃_ev_et_k)exp(-η̃_ev_eτ).

	
	Using the CFL condition max_e(v_e)τ/h≤1 leads to:
	
    [                                                                 S_1+Z_1                                                                       =         ∑_e=1^mv_e(∑_j=0^N-1(exp(-η_eh)-1)(f_e,j^k)^2p_eexp(-η_ex_e,j));                                                                                                                                               +                ∑_e=1^m(q_e^k)^2c_eexp(-η̃_ev_et_k)(exp(-η̃_ev_eτ)-1)1/τ;                                                                                                                                            CFL≤         ∑_e=1^mv_e(exp(-η_eh)-1)(∑_j=0^N-1(f_e,j^k)^2p_eexp(-η_ex_e,j));                                                                                                                                               +       ∑_e=1^m(q_e^k)^2c_eexp(-η̃_ev_et_k)(exp(-η̃_ev_eτ)-1)max_e(v_e)/h;                                                                                                                                               ≤    max_e(v_e(exp(-η_eh)-1))∑_e=1^m∑_j=0^N-1(f_e,j^k)^2p_eexp(-η_ex_e,j);                                                                                                                                               + max_e(v_e)max_e(exp(-η̃_ev_eτ)-1)∑_e=1^m(q_e^k)^2c_eexp(-η̃_ev_et_k)1/h;                                                                                                                                               ≤   1/hmax{max_e(v_e(exp(-η_eh)-1)),max_e(v_e)max_e(exp(-η̃_ev_eτ)-1)}V^k;                                                                                                                                               =                                                                 -ν V^k. ]

	Therefore, if the inequality assumption (<ref>) is fulfilled, i.e., S_2+Z_2≤ 0, we obtain:
	
    V^k+1-V^k/τ   ≤ S_1+Z_1
       ≤-ν V^k.

	Hence, we have V^k+1-V^k≤-τν V^k with ν=min{ν_1,ν_2}>0, whereby
	
    ν_1   =-1/hmax_e(v_e(exp(-η_eh)-1))=1/hmin_e(v_e(1-exp(-η_eh))) and
    ν_2   =-1/hmax_e(v_e)max_e(exp(-η̃_ev_eτ)-1)=1/hmax_e(v_e)min_e(1-exp(-η̃_ev_eτ)).

	Using the CFL condition max_e(v_e)τ/h≤1, implies
	
    [                                                                               0                                                                               <           τν=τ/hmin{min_e(v_e(1-exp(-η_eh))),max_e(v_e)min_e(1-exp(-η̃_ev_eτ))};                                                                                                                                                            CFL≤ 1/max_e(v_e)min{min_e(v_e(1-exp(-η_eh))),max_e(v_e)min_e(1-exp(-η̃_ev_eτ))}≤ 1. ]

	Finally, we prove estimate (<ref>) by recursively applying V^k+1≤(1-τν)V^k, which reads as
	
    V^k+1   ≤(1-τν)V^k≤(1-τν)^k+1V^0
       =exp((k+1)ln(1-τν))V^0≤exp(-ντ(k+1))V^0,

	where we use that ln(1-τν)≤-τν for τν∈[0,1].
	
	It remains to show (<ref>). It is well-known that the discrete Lyapunov function is bounded by the discrete norm in (<ref>):
	
    C̅(∑_j=0^N-1∑_e=1^m(f_e,j^k)^2h+∑_e=1^m(q_e^k)^2)≤ V^k≤C̃(∑_j=0^N-1∑_e=1^m(f_e,j^k)^2h+∑_e=1^m(q_e^k)^2),

	for some C̅,C̃>0.
	The above norm estimation in combination with (<ref>) leads to the desired result:
	
    ∑_j=0^N-1∑_e=1^m(f_e,j^k+1)^2h+∑_e=1^m(q_e^k+1)^2≤ Cexp(-ντ(k+1))(∑_j=0^N-1∑_e=1^m(f_e,j^0)^2h+∑_e=1^m(q_e^0)^2),

	for k∈{0,...,K-2} and C=C̃/C̅>0.


The inequality assumption (<ref>) in Theorem <ref> is a crucial point for exponential stability of the discretized problem. Therefore, we focus on this assumption in the following and start with the subsequent Remark <ref>.


	For sufficiently large k, assumption (<ref>) reduces to the condition that the matrix G^TP_0Λ G-P_N-1Λ is negative semidefinite, where G is the inflow matrix from (<ref>) and P_j is defined as P_j=diag(p_1exp(-η_1x_1,j),…,p_mexp(-η_mx_m,j)). 
	

	
	Let us consider again a serial system with two processing units and queues.
	The first intuition is then to insert the coupling conditions for f_e,-1^k. However, due to the queues determining the coupling conditions, this is 
	comparable to the non-queue case, cf. <cit.>. We can expect the queue at processor 2 to be equal to zero for t→ T, i.e., q_2^k=0 for k→(K-1), provided that T is large enough. Therefore, we differentiate two cases, i.e., a queue and no queue at processor 2. 
	
		
  * We obtain the non-queue case for k→(K-1) leading to the expression S_2+Z_2:
		
    [                                                                    S_2+Z_2                                                                     Z_2=0= ∑_e=1^2v_e((f_e,-1^k)^2p_eexp(-η_ex_e,0)-(f_e,N-1^k)^2p_eexp(-η_ex_e,N-1));                                                                                                                                                     =                     ∑_e=1^2v_e((∑_l=1^2G_e,lf_l,N-1^k)^2p_eexp(-η_ex_e,0));                                                                                                                                                     -                               ∑_e=1^2v_e((f_e,N-1^k)^2p_eexp(-η_ex_e,N-1));                                                                                                                                                     =                                  (f^k)^TG^TP_0Λ G(f^k)-(f^k)^TP_N-1Λ(f^k), ]

		with f^k=(f_1,N-1^k,f_2,N-1^k)^T. To satisfy assumption (<ref>), we have to ensure that the matrix G^TP_0Λ G-P_N-1Λ is negative semidefinite. Then, by assuming 
		
    G = [ 0 κ; 1 0 ],

		we are able to determine κ>0 such that a linear feedback control can be applied.
		
  * For small k we might have q_2^k>0 and f_2,-1^k=μ_2. Therefore, we are not able to find κ such that exponential stability is obtained. In fact, we have:
		
    [                                                                  S_2                                                                    = v_1((f_1,-1^k)^2p_1exp(-η_1 x_1,0)-(f_1,N-1^k)^2p_1exp(-η_1x_1,N-1));                                                                [1ex]                                                                    +     v_2((μ_2)^2p_2exp(-η_2 x_2,0)-(f_2,N-1^k)^2p_2exp(-η_2x_2,N-1)),;                                                          [1ex]
    		Z_2                                                             q_1^k=0=               2q_2^k(f_1,N-1^k-μ_2)c_2exp(-η̃_2v_2t_k)exp(-η̃_2v_2τ);                                                                [1ex]                                                                    +                 τ(f_1,N-1^k-μ_2)^2c_2exp(-η̃_2v_2t_k)exp(-η̃_2v_2τ). ]

		For simplicity, we assume η_e=η, η̃_e=η̃ and v_1=v_2=p_1=p_2=c_2=1.
		Hence, to ensure S_2+Z_2≤0, the feedback law u_1^k=f_1,-1^k from assumption (<ref>) has to fulfill:
		
    (u_1^k)^2   ≤(f_1,N-1^k)^2exp(-η x_1,N-1)+(f_2,N-1^k)^2exp(-η x_2,N-1)
       -(μ_2)^2exp(-η x_2,0)-2q_2^k(f_1,N-1^k-μ_2)exp(-η̃t_k)exp(-η̃τ)
       -τ(f_1,N-1^k-μ_2)^2exp(-η̃t_k)exp(-η̃τ) Y^k,

	
or, equivalently, if Y^k≥0, u_1^k≤√(Y^k).
	At the beginning, if queue 2 is not equal to zero yet, we use the time-dependent control law Y^k (<ref>). As soon as queue 2 is damped to 0, we can consider case 1 again. That means that for k→ (K-1), we are able to determine a matrix G, and thus a κ>0, such that exponential stability is ensured. 
	Consequently, it is not possible to find a general κ>0 for all t_k such that exponential stability is obtained. Quite the contrary, to include a linear feedback law, we have to differentiate between the two cases above. This explains the kink in the Lyapunov function when only the linear feedback (LF) with κ>0 for all time steps is used, see Figure <ref>. 
	
In the following, we call the feedback law that differentiates between the non-queue case and the case with queues the mixed feedback (MF).
We note that Y^k in (<ref>) 	coincides with the continuous version X(t) in (<ref>) (under the assumption that q_2^k>0 and v_e=1) as τ→0. 
A similar result is obtained for the numerical decay rate ν=min(ν_1,ν_2) (see (<ref>)) of the discrete Lyapunov function:

Assume τ = h/max_e(ṽ^e), i.e. CFL with equality, is satisfied. Then the numerical decay rate ν=min(ν_1,ν_2) converges to the analytical decay rate min(η,η̃)min_e(v_e) as h → 0
(see Lemma <ref>).


Since the minimum is continuous, we show the convergence of ν_1 →ηmin_e(v_e) and ν_2 →η̃min_e(v_e) separately.
We use a the Taylor expansion Tf(h;0)=η_e h+O(h^2) for f(h)1-exp(-η_eh) and obtain
	
    ν_1    = min_e(ṽ^e1/hTf(h;0)) = min_e(ṽ^e(η_e+O(h)))h→ 0⟶min_e(ṽ^eη_e).

Using the CFL condition with equality, we have 
	
    ν_2   =1/hmax_e(ṽ^e)min_e(1-exp(-η̃_eṽ^eτ)) =1/hmax_e(ṽ^e)min_e(1-exp(-η̃_eṽ^eh/max_e(ṽ^e))).

	We use the Taylor expansion Tf(h;0)=ṽ^e/max_e(ṽ^e)η̃_e h+O(h^2) for f(h)1-exp(-η̃_eṽ^eh/max_e(ṽ^e)). This yields
	
    ν_2    =min_e(1/hmax_e(ṽ^e)(1-exp(-η̃_eṽ^eh/max_e(ṽ^e))))
    
    	    = min_e(ṽ^e(η̃_e+O(h)))
    h→ 0⟶      min_e(ṽ^eη̃_e).


	




 §.§ Computational Experiments

We stick to the example presented in Figure <ref>.
First, we consider case 1, which means that we have to ensure that the matrix G^TP_0Λ G-P_N-1Λ is negative semidefinite. Since we consider a serial network, the matrix G is given by (<ref>).
The two eigenvalues of G^TP_0Λ G-P_N-1Λ are:

    λ_1   =p_2v_2exp(-η x_2,0)-p_1v_1exp(-η x_1,N-1),
    λ_2   =κ^2p_1v_1exp(-η x_1,0)-p_2v_2exp(-η x_2,N-1).

Assuming p_1=p_2=v_1=v_2=1, due to Theorem <ref>, we further have to consider the offset x_1,N-1=x_2,0=l to guarantee λ_1≤0:

    λ_1   =exp(-η l)-exp(-η l)=0,
    λ_2   =κ^2-exp(-η 2l).

To ensure the exponential stability, the eigenvalue λ_2 has to be nonpositive. Therefore, the linear feedback law is given by:

    u_1^k=κ f_2,N-1^k,
 
with κ>0 satisfying

    κ^2≤exp(-η2l)  ⇔ κ≤exp(-η l).
 

Note that the above case is only valid if the queue in front of processor 2 is equal to zero. If this is not the case, we implement the feedback law given by Y^k from (<ref>). Due to the offset in the spatial domain we obtain:

    [                                                                        (u_1^k)^2                                                                                ≤   (f_1,N-1^k)^2exp(-η x_1,N-1)+(f_2,N-1^k)^2exp(-η x_2,N-1)-(μ_2)^2exp(-η x_2,0);                                                                                                                                                                 - 2q_2^k(f_1,N-1^k-μ_2)exp(-η̃t_k)exp(-η̃τ)-τ(f_1,N-1^k-μ_2)^2exp(-η̃t_k)exp(-η̃τ);                                                                                                                                                           offset=                  (f_1,N-1^k)^2exp(-η l)+(f_2,N-1^k)^2exp(-η 2l)-(μ_2)^2exp(-η l);                                                                                                                                                                 - 2q_2^k(f_1,N-1^k-μ_2)exp(-η̃t_k)exp(-η̃τ)-τ(f_1,N-1^k-μ_2)^2exp(-η̃t_k)exp(-η̃τ);                                                                                                                                                                 =                                                                             Y^k. ]


So a possible choice for the mixed feedback law is given by:

    u_1^k=0,         if  Y^k<0,
    √(Y^k),    if Y^k≥ 0.

For the numerical computations, if not stated otherwise, we choose the parameter values η=η̃=1/2,l=1/2,h=0.01,T=30 and τ/hmax_ev_e=1. From (<ref>), we choose κ=e^-1/4≈0.7788. 
We set μ_1=6,μ_2=4 and choose the initial conditions f_1(0,x)=4,f_2(0,x)=4 and q_1(0)=0,q_2(0)=1. The velocities are v_e=1 for e∈{1,2}.
We denote the Lyapunov function, i.e. the upper bound in (<ref>), by:

    V_up^kexp(-ντ k)V^0.
 
In the first numerical example, we compare the linear feedback control from Section <ref>, with the mixed feedback control from Theorem <ref> so that there is a difference between the non-queue case and the case with positive queues.
Figure <ref> shows a decreasing Lyapunov function with a kink that results from using the linear feedback law with κ=0.5. 
Apparently, the discrete Lyapunov function lies underneath the upper bound V_up. However, the kink in the Lyapunov function can not be avoided. Due to the occurrence of the kink, the parameter value κ=0.7788 is not appropriate for the linear feedback since it would lead to a Lyapunov function that lies above V_up. For all values of κ, we can at most expect a decreasing discrete Lyapunov function, i.e., asymptotic stability.

For the mixed feedback law, the parameter value κ=0.7788 leads to the desired exponential decay of the Lyapunov function, see Figure <ref>, and V^k and V_up^k overlay almost exactly. The difference in the two approaches becomes even more clear if we have a look at the log-plot of both discrete Lyapunov functions. The log-plot of the discrete Lyapunov function of the mixed feedback is a straight line, whereas we can notice a kink in the log-plot for the linear feedback in Figure <ref>. 




For the same setting, we investigate the numerical behavior of the decay rate that is given by:

    ν=min{1/hmin_e(v_e(1-exp(-η_eh))),1/hmax_e(v_e)min_e(1-exp(-η̃_ev_eτ))}.

 To do so, we set η=η̃=0.575  and vary the value for the velocities. Table <ref> shows the convergence of the decay rate ν to min_e(v_e)min(η,η̃) for v_e=1 and v_e=1/2 as mentioned in Corollary <ref>. Furthermore, the table entries 
‖·‖_∞ and ‖·‖_L^2 denote the L^∞- and L^2-norm of the difference V_up^k-V^k. As expected, we observe first-order convergence of the discrete Lyapunov function towards the upper bound.

C[1]>m#1

The dependence of the discrete Lyapunov function on the parameter κ is shown in Table <ref>. Here, the ratio of the Lyapunov function at time t=T and t=0 is small for small values of κ and vice versa. The spatial step size h=0.00125 and v_e=v=1 imply that the decay rate ν is quite close to  η̃=min(η,η̃), see Table <ref>.


However, considering production systems with only decreasing queues does not seem reasonable. Therefore, we set μ_1=6, μ_2=4 and choose the initial conditions f_1(0,x)=6, f_2(0,x)=4 and q_2(0)=0, which lead to an increasing queue in front of processor 2. We set v_1=v_2=1, l=1/2, T=30 and η=η̃=0.2. The Lyapunov function V_1 for the flux and V_2 for the queue are depicted in Figure <ref>. 
We observe exponential decay of the composed Lyapunov function V=V_1+V_2, see Figure <ref>. Moreover, the outflow and the queue of processor 2 are damped to zero, see Figure <ref> and <ref>. The linear feedback takes over as soon as the queue is equal to zero. The reason for the shape of the feedback becomes more clear in the following example.



As a final experiment, we investigate the behavior of the Lyapunov function for a larger gap between the two maximal capacities μ_1 and μ_2. Therefore, we set μ_2=4 and vary the values for μ_1. We choose the initial conditions for the flux equal to the maximal capacities and assume again q_e(0)=0. 
As in the previous example, we set η=0.2, η̃=0.2, which leads to κ=e^-1/10≈0.9048. Regarding the linear feedback law, a higher μ_1 leads to a bigger kink in the Lyapunov function as shown in Figure <ref>. In comparison to that, we notice that the mixed feedback law drastically reduces the kink in the Lyapunov function and leads to the desired exponential decay. This shows that the mixed feedback law is even able to deal with higher queue loads. 

In Figure <ref>, the feedback u_1^k for both feedback laws is depicted. It should be noted that the mixed feedback u_1^k is decreasing more slowly as soon as the queue load is sinking (i.e. μ_2>f_1,N-1^k). The evolution of the queue is reflected in the feedback since it is part of the mixed feedback law, see (<ref>). This behavior becomes even more significant for larger μ_1 leading to a a higher queue load. 
Furthermore we remark that the mixed feedback slightly increases for all choices of μ_1 as soon as the linear feedback takes over due to the queue load being equal to zero. Shortly after the increase, the linear feedback starts to decrease again since the flux in the first processor has been already damped below μ_2.

Applying simply the linear feedback law leads to a constant feedback u_1^k=κ f_2,N-1^k=4κ at the beginning due to the queue being nonzero. 
The constant behavior of u_1^k remains maintained longer for larger μ_1, also resulting in a bigger kink of the Lyapunov function. Finally, both feedback laws damp the outflow and the queue to zero.







§ CONCLUSION

The main topic of this paper is the continuous and discrete feedback stabilization 
for a coupled PDE-ODE system describing a serial production system. 
By using of an adapted Lyapunov function, we are able to prove exponential stability 
and derive a linear feedback law. We also observe that the linear feedback law leads to at most asymptotic stability. 
Based on the numerical approximation for the coupled PDE-ODE system, we can compute decay rates
and derive a mixed feedback law which is consistent with the analytical result but able to resolve queueing situations. 




§ ACKNOWLEDGMENTS

This work was financially supported by BMBF project ENets (05M18VMA) and the DFG grants No. GO 1920/4-1 and GO 1920/7-1.



siam



