Feedback Stabilization for a coupled PDE-ODE Production System

TL;DR

This paper develops a Lyapunov-based control approach to stabilize a coupled PDE-ODE production system, ensuring exponential stability even in bottleneck scenarios through a mixed feedback law, supported by theoretical and numerical analysis.

Contribution

It introduces a novel feedback stabilization method for coupled PDE-ODE systems in production models, combining theoretical analysis with computational validation.

Findings

The mixed feedback law guarantees exponential stability.

Stability is achieved even in bottleneck conditions.

Numerical examples confirm theoretical results.

Abstract

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.

Tables

Table 1. Table 1 : Convergence of the decay rate ν 𝜈 \nu for η = η ~ = 0.575 𝜂 ~ 𝜂 0.575 \eta=\tilde{\eta}=0.575 and first-order convergence of the discretization for a CFL constant equal to 1,  N = 1 2 ​ h 𝑁 1 2 ℎ N=\frac{1}{2h} and different velocities

	$v_e=1$					$v_e=0.5$
$N$	${\parallel \cdot \parallel }_{\infty }$	Conv. Rate	${\parallel \cdot \parallel }_{L^2}$	Conv. Rate	$\nu$	${\parallel \cdot \parallel }_{\infty }$	Conv. Rate	${\parallel \cdot \parallel }_{L^2}$	Conv. Rate	$\nu$
10	0.0754	-	0.1326	-	0.5668	0.0834	-	0.2007	-	0.2834
50	0.0151	0.99	0.0265	1.00	0.5734	0.0153	1.05	0.0380	1.03	0.2867
100	0.0075	1.01	0.0132	1.00	0.5742	0.0076	1.01	0.0189	1.01	0.2871
200	0.0038	0.99	0.0066	1.00	0.5746	0.0038	1.00	0.0094	1.01	0.2873
400	0.0019	1.00	0.0033	1.00	0.5748	0.0019	1.00	0.0047	1.00	0.2874
800	0.0009	1.05	0.0017	0.97	0.5749	0.0009	1.08	0.0023	1.03	0.2874

Table 2. Table 2 : Dependence of the Lyapunov function on κ 𝜅 \kappa

$\kappa$	$\frac{V^T}{V^0}$	$\eta$	$\tilde{\eta }$	$\nu$
0.1	$3.75e^{-60}$	4.6052	0.5752	0.5750
0.25	$1.40e^{-36}$	2.7726	0.5752	0.5750
0.5	$1.05e^{-18}$	1.3863	0.5752	0.5750
0.75	$3.20e^{-8}$	0.5752	0.5752	0.5750

Equations

$${\partial }_t{\rho }_e(t,x)+{\partial }_{x}f({\rho }_e(t,x))=0,$$

$$f_e({\rho }_e(t,x))=\min \{v_e{\rho }_e(t,x),{\mu }_e\},\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\forall x\in [a_e,b_e],t\in {\mathbb{R}}_{+},$$

$${\rho }_e(0,x)={\rho }_{e,0}(x),$$

$$g_{in,e}(t):=f_{e-1}({\rho }_{e-1}(t,b_{e-1})),\text{leavevmode}\text{forleavevmodenobreak }e=2,...,m.$$

$${\partial }_t{q}_e(t)=g_{in,e}(t)-g_{out,e}(t),\text{leavevmode}\text{withleavevmodenobreak }{g}_{out,e}(t):=f_e({\rho }_e(t,a_e))$$

$$g_{out,e}(t):=f_e({\rho }_e(t,a_e))=\{\begin{array}{ll}\min \{f_{e-1}({\rho }_{e-1}(t,b_{e-1})),{\mu }_e\},& q_e(t)=0,\\ {\mu }_e,& q_e(t)>0.\end{array}$$

$$f_e({\rho }_e(t,a_e))\approx \min \left\{{\mu }_e,\frac{q_e(t)}{ϵ}\right\},\text{leavevmode}\text{with}\text{leavevmode}ϵ\ll 1.$$

$${\partial }_t{\rho }_e+v_e{\partial }_x{\rho }_e=0.$$

$$\Lambda =\text{diag}(v_1,\dots ,v_m),f(t,x)=(f_1(t,x),\dots ,f_m(t,x){)}^T,{\partial }_{t}q(t)=({\partial }_t{q}_1(t),\dots ,{\partial }_t{q}_m(t){)}^T,$$

$$g_{in}(t)=(g_{in,1}(t),\dots ,g_{in,m}(t){)}^T\text{ and }{g}_{out}(t)=(g_{out,1}(t),\dots ,g_{out,m}(t){)}^T,$$

$${\partial }_{t}f(t,x)+\Lambda {\partial }_{x}f(t,x)=0,\text{leavevmode}x\in [0,l],t\in {\mathbb{R}}_{+},$$

$${\partial }_{t}q(t)=g_{in}(t)-g_{out}(t),$$

$$\bar{v}\frac{\tau }{h}\le 1,\text{leavevmode}\text{leavevmode}\text{leavevmode}\bar{v}:=\underset{e\in \{1,...,m\}}{\max }{v}_e.$$

$$f_{e,j}^{k+1}$$

$$f_{e,-1}^k$$

$$f_{1,-1}^k$$

$$q_e^{k+1}$$

$${g_{in,e}}^k$$

$$f_{e,j}^0$$

$$\Vert (f(t,\cdot ),q(t))\Vert =\sqrt{\Vert f(t,\cdot ){\Vert }_{L^2((0,l);{\mathbb{R}}^m)}^2+|q(t){|}^2},$$

$$\Vert f(t,\cdot ){\Vert }_{L^2((0,l);{\mathbb{R}}^m)}^2+|q(t){|}^2\le C{e}^{-\nu t}\left(\Vert f(0,\cdot ){\Vert }_{L^2((0,l);{\mathbb{R}}^m)}^2+|q(0){|}^2\right),\text{leavevmode}\forall t\in {\mathbb{R}}_{+}.$$

$$z=\left(z_1,\dots ,z_m\right),\text{leavevmode}\text{with}\text{leavevmode}{z}_e=\zeta \underset{e}{\min }({\mu }_e)\text{leavevmode}\text{for}\text{leavevmode}0\le \zeta <1.$$

$$V(t)$$

$$\displaystyle=\int_{0}^{l}f^{T}(t,x)P(x)f(t,x)\leavevmode\nobreak\ dx+q^{T}(t)Q(t)q(t),\leavevmode\nobreak\$$

$$V^k=\sum _{j=0}^{N-1}\sum _{e=1}^m(f_{e,j}^k{)}^2{p}_e{e}^{-{\eta }_e{x}_{e,j}}h+\sum _{e=1}^m(q_e^k{)}^2{c}_e{e}^{-{\tilde{\eta }}_e{v}_e{t}_k},\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{with}\text{leavevmode}{\eta }_e,{\tilde{\eta }}_e,p_e,c_e>0.$$

$$f_{1,-1}^k$$

$$q_e^{k+1}$$

$${g_{in,e}}^k$$

$$\dot{V}\le -\nu V,$$

$$-f(t,l{)}^T\Lambda P(l)f(t,l)+f(t,0{)}^T\Lambda P(0)f(t,0)+2{\partial }_{t}q(t{)}^{T}Q(t)q(t)\le 0$$

$$\nu =\min (\eta ,\tilde{\eta })\underset{e}{\min }(v_e)>0.$$

$$\begin{array}[b]{r@{}c@{}l}\dot{V}&\overset{\eqref{eq:lyap}}{=}&\int_{0}^{l}\left(\partial_{t}f^{T}P(x)f+f^{T}P(x)\partial_{t}f\right)dx+2\partial_{t}q^{T}Q(t)q+q^{T}\partial_{t}Q(t)q\\[4.30554pt] &\overset{\eqref{eq:fcompact1}}{=}&\int_{0}^{l}\left(-\partial_{x}f^{T}\Lambda P(x)f-f^{T}P(x)\Lambda\partial_{x}f\right)dx+2\partial_{t}q^{T}Q(t)q-\tilde{\eta}q^{T}\Lambda Q(t)q\\[4.30554pt] &=&-\int_{0}^{l}\partial_{x}\left[f^{T}\Lambda P(x)f\right]dx-\eta\int_{0}^{l}f^{T}\Lambda P(x)fdx+2\partial_{t}q^{T}Q(t)q-\tilde{\eta}q^{T}\Lambda Q(t)q\\[4.30554pt] &\overset{}{\leq}&-\left[f^{T}\Lambda P(x)f\right]_{0}^{l}-\eta\min_{e}({v}_{e})\int_{0}^{l}f^{T}P(x)fdx+2\partial_{t}q^{T}Q(t)q-\tilde{\eta}\min_{e}({v}_{e})q^{T}Q(t)q\\[4.30554pt] &\leq&-\min(\eta,\tilde{\eta})\min_{e}({v}_{e})V(t)-\left[f^{T}\Lambda P(x)f\right]_{0}^{l}+2\partial_{t}q^{T}Q(t)q.\end{array}$$

$$-{\left[f^T\Lambda P(x)f\right]}_0^l+2{\partial }_t{q}^{T}Q(t)q\le 0.\text{qed}$$

$$-v_1{f}_1(t,l{)}^2{p}_1{e}^{-\eta l}+v_1{f}_1(t,0{)}^2{p}_1{e}^{-\eta 0}-v_2{f}_2(t,l{)}^2{p}_2{e}^{-\eta l}+v_2{f}_2(t,0{)}^2{p}_2{e}^{-\eta 0}$$

$$+2q_1(t)(g_{in,1}(t)-g_{out,1}(t))c_1{e}^{-\tilde{\eta }{v}_{1}t}+2q_2(t)(g_{in,2}(t)-g_{out,2}(t))c_2{e}^{-\tilde{\eta }{v}_{2}t}\le 0.$$

$$v_1{u}_1(t{)}^2{p}_1$$

$$-v_2{f}_2(t,0{)}^2{p}_2-2q_2(t)g_{in,2}(t)c_2{e}^{-\tilde{\eta }{v}_{2}t}+2q_2(t)g_{out,2}(t)c_2{e}^{-\tilde{\eta }{v}_{2}t}.$$

Full text

AMS Classification: 65Mxx, 93D05, 90B30

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

1 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 [2, 3, 9, 14, 15, 16, 17]. The extension to networks has been originally introduced in [19], 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 [16] 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 [5, 10, 21, 23, 25] has been investigated intensively. However, there is only little literature available on the stability analysis for nonlinear coupled systems, see for example [24]. 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 [6, 11, 22]. 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 [12, 13]. Also, explicit decay rates for numerical schemes have only recently been established. In [4], exponential decay on a finite time horizon has been developed for non-conservative (first-order) schemes and in [20] 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 2, we introduce the production model that presented in [16] and provide a suitable numerical discretization. We make use of the theory of Lyapunov functions (Section 3.1) to derive a feedback law that ensures the exponential stability for the coupled PDE-ODE model in Section 3.2. 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 4. On the discrete level, we obtain a so-called mixed feedback law to deal with queueing situations, i.e. non-empty queue loads.

2 Model Equations and Numerical Discretization

First, we recall the production model taken from [16]. 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\in\{1,...,m\}$. Furthermore, the maximal capacity $\mu_{e}>0$ and the processing velocity $v_{e}>0$ of each unit are constant parameters. The dynamics of the product density $\rho_{e}(t,x)$ defined on a segment $[a_{e},b_{e}]$ follows [2]:

[TABLE]

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:

[TABLE]

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 $\mu_{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:

[TABLE]

describing the outflow from queue $e$ to processor $e$ at time $t$, cf. Figure 1.

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 $\mu_{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

[TABLE]

Equations (2.3) and (2.4) ensure that the production model is well-defined, provided that $f_{in}(t)$ and $\rho_{e,0}(x)$ are of total bounded variation, see [16].

To avoid discontinuities in the derivative of the queue length, a smoothed out version of (2.4) has been derived in [1] and reads as

[TABLE]

Provided ${v}_{e}\rho_{e,0}\leq\mu_{e}$, (2.1a) reduces to a linear advection equation

[TABLE]

For the numerical investigations later on, we introduce a first-order discretization of the coupled production model. From now on, we assume that $x\in[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 3.1. By defining

[TABLE]

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

[TABLE]

equipped with initial values $f_{e,0}(x)=v_{e}\rho_{e,0}(x),\leavevmode\nobreak\ q(0)=q_{0}\geq 0$. Note that (2.6a) is already given in characteristic form, see [6].

For the discretization of (2.6a) 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 [16]. We use an equidistant grid with constant time step size $\tau$ 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+\frac{1}{2})h$ for $j\in\{-1,\ldots,N-1\}$. The interfaces of the cells are located at $x_{j-\frac{1}{2}}=jh$ for $j\in\{-1,...,N-1\}$. The considered time domain is $[0,T]$, whereby the discrete time is denoted by $t_{k}=k\tau$ for $k\in\{0,...,K-1\}$ such that $(K-1)\tau=T$. Moreover, the CFL condition is assumed to be fulfilled, i.e.,

[TABLE]

For further investigations, we assume that $q_{1}^{k}=0\leavevmode\nobreak\ \forall k$, which is satisfied if $f_{in}^{k}\leq\mu_{1}\leavevmode\nobreak\ \forall k$. Then, the discretized equations for $j\in\{0,...,N-1\}$ and $k\in\{0,...,K-2\}$ read:

[TABLE]

Here, (2.8b),(2.8c) 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,\dots,m$). Due to assumption $f_{in}^{k}\leq\mu_{1}\leavevmode\nobreak\ \forall 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.

3 Feedback Stabilization for the Continuous System

We now investigate the stabilization of the continuous model (2.6b). 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.

3.1 Lyapunov Stability

In the following, we analyze the exponential stability of the continuous system (2.6b). In contrast to [6, 8], we are faced with queues and therefore investigate the exponential stability in the norm

[TABLE]

where $|\cdot|$ is the Euclidean norm. The introduced norm (3.1) is the natural choice of the product space $L^{2}((0,l);\mathbb{R}^{m})\times\mathbb{R}^{m}$ motivated by [21, 23].

Definition 3.1**.**

The system (2.6a)-(2.6b) is exponentially stable in the sense of the norm (3.1), if there exist $\nu>0$ and $C>0$ such that for every initial condition $f_{0}\in L^{2}((0,l);\mathbb{R}^{m})$ and $q(0)\in\mathbb{R}^{m}$ the solution to (2.6b) satisfies

[TABLE]

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

[TABLE]

By translation $\tilde{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. [6, 18, 20]) needs to be adapted. We define the following Lyapunov function candidate:

[TABLE]

with the weighting matrices $P(x)=\text{diag}(p_{1}e^{-\eta_{1}x},...,p_{m}e^{-\eta_{m}x})$, $p_{e}>0,\eta_{e}>0$ and

$Q(t)=\text{diag}(c_{1}e^{-\tilde{\eta}_{1}{v}_{1}t},\ldots,c_{m}e^{-\tilde{\eta}_{m}{v}_{m}t})$, $c_{e}>0,\tilde{\eta}_{e}>0$. In order to include the ODEs for the queues, we add the quadratic form $V_{2}(t)=q(t)^{T}Q(t)q(t)$. If not stated otherwise, we assume $\eta_{e}=\eta$ and $\tilde{\eta}_{e}=\tilde{\eta}$, independent of the individual properties.

Remark 3.1**.**

For the discrete stability analysis in Section 4, we define the discrete version of the Lyapunov function in (3.3) as

[TABLE]

3.2 Analytical Feedback Stabilization

The goal is to derive a feedback law using the Lyapunov function from (3.3). Figure 2 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.

Remark 3.2**.**

By considering a control input $u_{1}(t)$ as an inflow into the system, equations (2.8c) and (2.8d) of the numerical discretization are adapted as follows:

[TABLE]

Since so far, only the non-queue case has been studied in the literature [6, 7, 8]. 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 (3.3).

Lemma 3.1**.**

The continuous solution satisfies

[TABLE]

i.e., exponential decay of $V$ is obtained, if

[TABLE]

is fulfilled $\forall t\in\mathbb{R_{+}}$. The decay rate $\nu$ is given by

[TABLE]

Proof.

In the following, we assume that $f$ and $q$ are sufficiently smooth. We obtain:

[TABLE]

To obtain the exponential decay of $V$, we deduce the condition

[TABLE]

Next, we derive a feedback law by using Lemma 3.7 and by imposing the feedback control $u_{1}(t)$ at processor 1. For the serial network in Figure 2 we have to show

[TABLE]

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

[TABLE]

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 (2.5) for $\epsilon\ll 1$ to ensure enough regularity:

[TABLE]

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\equiv 0$ and $q\equiv 0$, the control $u_{1}=0$ is obtained. Assuming ${v}_{e}={v}$ for $e\in\{1,2\}$, we get:

[TABLE]

Assuming ${v}_{e}={v}$ and $p_{e}=p=c_{e}=c=1$, this result can be extended to a serial system with $m$ processors:

[TABLE]

If $X(t)<0$, we may set $u_{1}(t)=0$. Otherwise, we set $u_{1}(t)=\sqrt{X(t)}$ motivated by (3.9). 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,

[TABLE]

for $\kappa>0$.

In the case of vanishing queues, i.e. $q_{e}(t)=0$ for all $t\in\mathbb{R}_{+}$, the results from Lemma 3.7 and (3.9)–(3.12) coincide with [6, 8, 20] by a direct comparison.

3.3 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=\tau=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 (3.12) 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 $\mu_{1}=10,\leavevmode\nobreak\ \mu_{2}=9,\leavevmode\nobreak\ \mu_{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 $\eta_{e}=\tilde{\eta}_{e}=0.1$ and $p_{e}=c_{e}=1$. We denote the upper bound from the proof of Lemma 3.7 by:

[TABLE]

Figure 3(a) shows that the estimated Lyapunov function for $\kappa=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 3(b).

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 4.

4 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 [4, 20]. 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 (3.4) also requires the definition of a discrete norm. Therefore, we use for the stability analysis of the disctrized equations

[TABLE]

where $\lVert f\rVert^{2}_{h,2}=\sum\limits_{j=0}^{N-1}\sum\limits_{e=1}^{m}(f_{e,j})^{2}h$ is the discrete $L^{2}$-norm of $f$.

Since we are interested in a feedback stabilization for the numerical discretization, we equip (2.6b) with

[TABLE]

where the matrix $G\in\mathbb{R}^{m\times 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.

Theorem 4.1**.**

Let $T>0$ be arbitrarily large but fixed, $f_{in}^{k}\leq\mu_{1}$ and the feedback law $f_{in}^{k}=u_{1}^{k}=f_{1,-1}^{k}$ satisfies for $k\in\{0,...,K-1\}$ the inequality:

[TABLE]

Then, provided the CFL condition (2.7) is fulfilled, the following holds: there exist $\eta_{e}>0$ and $\tilde{\eta}_{e}>0$ such that the numerical scheme (2.8a)-(2.8e) is exponentially stable, i.e., the numerical solution satisfies

[TABLE]

where the decay rate $\nu$ is given by:

[TABLE]

Furthermore, $f_{e,j}^{k+1}$ and $q_{e}^{k+1}$ are exponentially stable in the sense of the discrete norm in (4.1), i.e.,

[TABLE]

for $k\in\{0,...,K-2\}$ and $C>0$.

Proof.

In order to obtain the exponential stability of the numerical solution, we show that the discrete time derivative of the discrete Lyapunov function in (3.4) fulfills:

[TABLE]

By inserting the discrete Lyapunov function at time $t^{k}$ and $t^{k+1}$, we obtain:

[TABLE]

with

[TABLE]

At first, we have a closer look at $C_{1}$ by using the numerical discretization scheme (2.8a)-(2.8c) for the advection equation:

[TABLE]

Applying the CFL condition, we obtain

[TABLE]

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:

[TABLE]

Therefore, we have

[TABLE]

with

[TABLE]

Regarding $C_{2}$, using the numerical discretization, or more specifically, the explicit Euler method for the queues, we obtain:

[TABLE]

In the next step, we use the binomial formula, which yields:

[TABLE]

with

[TABLE]

Using the CFL condition $\max_{e}({v}_{e})\frac{\tau}{h}\leq 1$ leads to:

[TABLE]

Therefore, if the inequality assumption (4.3) is fulfilled, i.e., $S_{2}+Z_{2}\leq 0$, we obtain:

[TABLE]

Hence, we have $V^{k+1}-V^{k}\leq-\tau\nu V^{k}$ with $\nu=\min\{\nu_{1},\nu_{2}\}>0$, whereby

[TABLE]

Using the CFL condition $\max_{e}({v}_{e})\frac{\tau}{h}\leq 1$, implies

[TABLE]

Finally, we prove estimate (4.4) by recursively applying $V^{k+1}\leq(1-\tau\nu)V^{k}$, which reads as

[TABLE]

where we use that $\ln(1-\tau\nu)\leq-\tau\nu$ for $\tau\nu\in[0,1]$.

It remains to show (4.6). It is well-known that the discrete Lyapunov function is bounded by the discrete norm in (4.1):

[TABLE]

for some $\bar{C},\tilde{C}>0$. The above norm estimation in combination with (4.9) leads to the desired result:

[TABLE]

for $k\in\{0,...,K-2\}$ and $C=\frac{\tilde{C}}{\bar{C}}>0$. ∎

The inequality assumption (4.3) in Theorem 4.1 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 4.1.

Remark 4.1**.**

For sufficiently large $k$, assumption (4.3) reduces to the condition that the matrix $G^{T}P_{0}\Lambda G-P_{N-1}\Lambda$ is negative semidefinite, where $G$ is the inflow matrix from (4.2a) and $P_{j}$ is defined as $P_{j}=\text{diag}(p_{1}\exp(-\eta_{1}x_{1,j}),\ldots,p_{m}\exp(-\eta_{m}x_{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. [20]. We can expect the queue at processor 2 to be equal to zero for $t\rightarrow T$, i.e., $q_{2}^{k}=0$ for $k\rightarrow(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\rightarrow(K-1)$ leading to the expression $S_{2}+Z_{2}$:

[TABLE]

with $f^{k}=(f_{1,N-1}^{k},f_{2,N-1}^{k})^{T}$. To satisfy assumption (4.3), we have to ensure that the matrix $G^{T}P_{0}\Lambda G-P_{N-1}\Lambda$ is negative semidefinite. Then, by assuming

[TABLE]

we are able to determine $\kappa>0$ such that a linear feedback control can be applied. 2. 2.

For small $k$ we might have $q_{2}^{k}>0$ and $f_{2,-1}^{k}=\mu_{2}$. Therefore, we are not able to find $\kappa$ such that exponential stability is obtained. In fact, we have:

[TABLE]

For simplicity, we assume $\eta_{e}=\eta$, $\tilde{\eta}_{e}=\tilde{\eta}$ and ${v}_{1}={v}_{2}=p_{1}=p_{2}=c_{2}=1$. Hence, to ensure $S_{2}+Z_{2}\leq 0$, the feedback law $u_{1}^{k}=f_{1,-1}^{k}$ from assumption (4.3) has to fulfill:

[TABLE]

or, equivalently, if $Y^{k}\geq 0$, $u_{1}^{k}\leq\sqrt{Y^{k}}$. At the beginning, if queue 2 is not equal to zero yet, we use the time-dependent control law $Y^{k}$ (4.11). As soon as queue 2 is damped to 0, we can consider case 1 again. That means that for $k\rightarrow(K-1)$, we are able to determine a matrix $G$, and thus a $\kappa>0$, such that exponential stability is ensured. Consequently, it is not possible to find a general $\kappa>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 $\kappa>0$ for all time steps is used, see Figure 3(a).

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 $\eqref{eq:u1k}$ coincides with the continuous version $X(t)$ in (3.10) (under the assumption that $q_{2}^{k}>0$ and ${v}_{e}=1$) as $\tau\rightarrow 0$. A similar result is obtained for the numerical decay rate $\nu=\min(\nu_{1},\nu_{2})$ (see (4.5)) of the discrete Lyapunov function:

Corollary 4.1**.**

Assume $\tau=\frac{h}{\max_{e}(\tilde{v}^{e})}$, i.e. CFL with equality, is satisfied. Then the numerical decay rate $\nu=\min(\nu_{1},\nu_{2})$ converges to the analytical decay rate $\min(\eta,\tilde{\eta})\min_{e}({v}_{e})$ as $h\to 0$ (see Lemma 3.7).

Proof.

Since the minimum is continuous, we show the convergence of $\nu_{1}\to\eta\min_{e}({v}_{e})$ and $\nu_{2}\to\tilde{\eta}\min_{e}({v}_{e})$ separately. We use a the Taylor expansion $Tf(h;0)=\eta_{e}h+O(h^{2})$ for $f(h)\coloneqq 1-\exp(-\eta_{e}h)$ and obtain

[TABLE]

Using the CFL condition with equality, we have

[TABLE]

We use the Taylor expansion $Tf(h;0)=\frac{\tilde{v}^{e}}{\max_{e}(\tilde{v}^{e})}\tilde{\eta}_{e}h+O(h^{2})$ for $f(h)\coloneqq 1-\exp\left(-\tilde{\eta}_{e}\tilde{v}^{e}\frac{h}{\max_{e}(\tilde{v}^{e})}\right)$. This yields

[TABLE]

∎

4.1 Computational Experiments

We stick to the example presented in Figure 1. First, we consider case 1, which means that we have to ensure that the matrix $G^{T}P_{0}\Lambda G-P_{N-1}\Lambda$ is negative semidefinite. Since we consider a serial network, the matrix $G$ is given by (4.10). The two eigenvalues of $G^{T}P_{0}\Lambda G-P_{N-1}\Lambda$ are:

[TABLE]

Assuming $p_{1}=p_{2}={v}_{1}={v}_{2}=1$, due to Theorem 4.1, we further have to consider the offset $x_{1,N-1}=x_{2,0}=l$ to guarantee $\lambda_{1}\leq 0$:

[TABLE]

To ensure the exponential stability, the eigenvalue $\lambda_{2}$ has to be nonpositive. Therefore, the linear feedback law is given by:

[TABLE]

with $\kappa>0$ satisfying

[TABLE]

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 (4.11). Due to the offset in the spatial domain we obtain:

[TABLE]

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

[TABLE]

For the numerical computations, if not stated otherwise, we choose the parameter values $\eta=\tilde{\eta}=\frac{1}{2},l=\frac{1}{2},h=0.01,T=30$ and $\frac{\tau}{h}\max_{e}{v}_{e}=1$. From (4.12), we choose $\kappa=e^{-\frac{1}{4}}\approx 0.7788$. We set $\mu_{1}=6,\mu_{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\in\{1,2\}$. We denote the Lyapunov function, i.e. the upper bound in (4.4), by:

[TABLE]

In the first numerical example, we compare the linear feedback control from Section 3.2, with the mixed feedback control from Theorem 4.1 so that there is a difference between the non-queue case and the case with positive queues. Figure 4(a) shows a decreasing Lyapunov function with a kink that results from using the linear feedback law with $\kappa=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 $\kappa=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 $\kappa$, we can at most expect a decreasing discrete Lyapunov function, i.e., asymptotic stability.

For the mixed feedback law, the parameter value $\kappa=0.7788$ leads to the desired exponential decay of the Lyapunov function, see Figure 4(b), 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 5(a).

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

[TABLE]

To do so, we set $\eta=\tilde{\eta}=0.575$ and vary the value for the velocities. Table 1 shows the convergence of the decay rate $\nu$ to $\min_{e}({v}_{e})\min(\eta,\tilde{\eta})$ for ${v}_{e}=1$ and ${v}_{e}=\frac{1}{2}$ as mentioned in Corollary 4.1. Furthermore, the table entries $\lVert\cdot\rVert_{\infty}$ and $\lVert\cdot\rVert_{L^{2}}$ denote the $L^{\infty}$- 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.

The dependence of the discrete Lyapunov function on the parameter $\kappa$ is shown in Table 2. Here, the ratio of the Lyapunov function at time $t=T$ and $t=0$ is small for small values of $\kappa$ and vice versa. The spatial step size $h=0.00125$ and ${v}_{e}={v}=1$ imply that the decay rate $\nu$ is quite close to $\tilde{\eta}=\min(\eta,\tilde{\eta})$, see Table 2.

However, considering production systems with only decreasing queues does not seem reasonable. Therefore, we set $\mu_{1}=6$, $\mu_{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=\frac{1}{2}$, $T=30$ and $\eta=\tilde{\eta}=0.2$. The Lyapunov function $V_{1}$ for the flux and $V_{2}$ for the queue are depicted in Figure 6(b). We observe exponential decay of the composed Lyapunov function $V=V_{1}+V_{2}$, see Figure 6(a). Moreover, the outflow and the queue of processor 2 are damped to zero, see Figure 6(c) and 6(d). 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 $\mu_{1}$ and $\mu_{2}$. Therefore, we set $\mu_{2}=4$ and vary the values for $\mu_{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 $\eta=0.2$, $\tilde{\eta}=0.2$, which leads to $\kappa=e^{-\frac{1}{10}}\approx 0.9048$. Regarding the linear feedback law, a higher $\mu_{1}$ leads to a bigger kink in the Lyapunov function as shown in Figure 7(a). 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 7(b), 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. $\mu_{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 (4.11). This behavior becomes even more significant for larger $\mu_{1}$ leading to a a higher queue load. Furthermore we remark that the mixed feedback slightly increases for all choices of $\mu_{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 $\mu_{2}$.

Applying simply the linear feedback law leads to a constant feedback $u_{1}^{k}=\kappa f_{2,N-1}^{k}=4\kappa$ at the beginning due to the queue being nonzero. The constant behavior of $u_{1}^{k}$ remains maintained longer for larger $\mu_{1}$, also resulting in a bigger kink of the Lyapunov function. Finally, both feedback laws damp the outflow and the queue to zero.

5 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.