Optimal control of a 2D diffusion-advection process with a team of mobile actuators under jointly optimal guidance^footnoteinfo

Recent development of mobile robots (unmanned aerial vehicles, terrestrial robots, and underwater vehicles) has greatly extended the type of distributed parameter system (DPS) over which mobile actuation and sensing can be deployed. Such a system is often modeled by a partial differential equation (PDE), which varies in both time and space. Exemplary applications of mobile control and estimation of a DPS can be found in thermal manufacturing [13], monitoring and neutralizing groundwater contamination [12], and wildfire monitoring [19].

We propose an optimization framework that simultaneously solves for the guidance of a team of mobile actuators and the control of a DPS. We consider a 2D diffusion-advection process as the DPS for its capability of modeling a variety of processes governed by continuum mechanics and the convenience of the state-space representation. The framework minimizes an integrated cost function, evaluating both the control of the DPS and the guidance of the actuators, subject to the dynamics of the DPS and the mobile actuators. The problem addresses the mobile actuator and the DPS as a unified system, instead of solely controlling the DPS. Furthermore, the additional degree of freedom endowed by mobility yields improved control performance in comparison to using stationary actuators.

The study of the control of a PDE-modeled DPS dates back to the 1960s (see the survey [29]). For fundamental results, see the textbooks [7, 28, 1]. Although it is possible to categorize prior work by whether the input operator is bounded, our literature review is categorized by the location of the actuation. When actuation acts on the boundary of the spatial domain, it is called boundary control. The main complexity in boundary control is that the input operator, which associates with the PDE’s boundary condition, is unbounded. This is addressed in the tutorial [15]. Recent developments on the design of boundary control uses the backstepping method [33], where a Volterra transformation determines a stabilizing control by transforming the system into a stable target system. When actuation acts in the interior of the spatial domain, it is called distributed control. For distributed control, the DPS is actuated by in-domain actuators that are either stationary or mobile.

The problem of determining the location of stationary actuators is called the actuator placement problem. Actuator placement has been studied for optimality in the sense of linear-quadratic (LQ) [24], H₂ [25], and H_∞ [17]. The author of [24] studies the actuator placement problem with the LQ performance criterion. The actuators’ locations are chosen to minimize the operator norm or trace norm of the Riccati operator solved from an algebraic Riccati equation associated with an infinite-dimensional system. If the input operator is compact and continuous with respect to actuator location in the operator norm, then a solution exists for the problem minimizing the operator norm of the Riccati operator [24, Theorem 2.6], under stabilizability and detectability assumptions. When computing the optimal actuator locations, if the approximated Riccati operator converges to the original Riccati operator at each actuator location, then the approximate optimal actuator locations converge to the exact optimal locations [24, Theorem 3.5]. For the above results to hold when minimizing the trace-norm of the Riccati operator, the input and output spaces have to be finite-dimensional [24, Theorems 2.10 and 3.9] in addition to the assumptions stated above.

The authors of [25] design optimal actuator placement by minimizing the H₂-control performance criterion, which minimizes the L₂-norm of the linear output of a linear system, subject state disturbances. Roughly speaking, H₂-control performance reduces the response to the disturbances while setting a zero initial condition, whereas the LQ performance reduces the response to the initial condition in a disturbance-free setting. For disturbances with known or unknown spatial distribution, the trace of the Riccati solution (scaled by the disturbance’s spatial distribution) or operator norm of the Riccati solution are minimized, respectively, where the existence of a solution and convergence to the exact optimal solution of the approximate optimal solution are guaranteed. In [17], the H_∞-performance criterion is minimized for actuator placement. Specifically, the actuators are placed in the locations that yield infimum of the attenuation bound (upper bound of the infinity norm of the closed-loop disturbance-to-output transfer function). The conditions for the existence of a solution and convergence to the exact optimal placement of the approximate optimal placement are provided.

Geometric approaches have also been proposed for actuator placement. For example, a modified centroidal Voronoi tessellation (mCVT) yields locations of actuators and sensors that yields least-squares approximate control and estimation kernels for a parabolic PDE [10]. The input operator is designed by maximizing the H₂-norm of the input-to-state transfer function, whereas the kernel of the state feedback is obtained using the Riccati solution. Next, mCVT determines the partition such that the actuator and sensor locations achieve optimal performance (in the sense of least-squares) to the input operator and state feedback kernel, respectively. A comparison of various performance criteria for actuator placement has been conducted for controlling a simply supported beam [26] and a diffusion process [27]. It has been analyzed that maximizing the minimum eigenvalue of the controllability gramian is not a useful criterion. Because the lower bound of the eigenvalue is zero, the minimum eigenvalue approaches zero as the dimension of approximation increases [26, 27].

The guidance of mobile actuators is designed to improve the control performance in comparison to stationary actuators. Various performance criteria have been proposed for guidance. In [13], a mobile heat source is steered to maintain a spatially uniform temperature distribution in 1D using the optimal control method. The formulation uses a finite-dimensional approximation for modeling the process and evaluating performance. Additionally, the admissible locations of the heat source are chosen within a discrete set that yields approximate controllability requirements. Algorithms are provided to solve the proposed problem with considerations on real-time implementation and hardware constraints. Experimental results demonstrate the performance of the proposed scheme. The authors of [14] propose an optimization framework that steers mobile actuators to control a reaction-diffusion process. A specific cost function, consisting of norms of control input and measurement of the DPS and the steering force, is minimized subject to dynamics of the actuator’s motion and the PDE, and bounds on the control input and state of the DPS. The implementation of the framework is emphasized by discrete mechanics and model predictive control which yield computationally tractable solutions, in addition to an approximation of the PDE and a discrete set of admissible actuator locations.

The problem of ultraviolet curing using a mobile radiant actuator is investigated in [38], where the curing process is modeled by a 1D nonlinear PDE. Both the radiant power and scanning velocity of the actuator are computed for reaching a target curing state. A dual extended Kalman filter is applied to estimate the state and parameters of the curing process for feedback control, based on the phases of curing. In [11], a navigation problem is studied in which a mobile agent moves through a diffusion process represented by a hazardous field with given initial and terminal positions. Such a formulation may be applied to emergency evacuation guidance from an indoor environment with carbon monoxide. Both problems with minimum time and minimum accumulated effects of hazards are formulated, and closed-form solutions are derived using the Hamiltonian. A Lyapunov-based guidance strategy for collocated sensors and actuators to estimate and control a diffusion process is proposed in [8]. The decentralized guidance of the actuators for controlling a diffusion process to a zero state is derived by constructing suitable Lyapunov functions. The same methodology is applied to construct a distributed consensus filter via the network among agents to improve state estimation. A follow-up work [9] incorporates nonlinear vehicle dynamics in addition.

The problem formulation in this paper includes a cost function that simultaneously evaluates controlling the PDE-modeled DPS, referred to as the PDE cost, and steering the mobile actuators, referred to as the mobility cost. The PDE cost is a quadratic cost of the PDE state and control, whose optimal value can be obtained by solving an operator-valued differential Riccati equation. Our results are based on the related work [2], which establishes Bochner integrable solutions of finite-horizon Riccati integral equations (with values in Schatten p-classes) associated with infinite-dimensional systems. The existence conditions for the solution of exact and approximate Riccati integral equations are established in [2]. The significance of the Bochner integrable solution is that it allows the implementation of simple numerical quadratures for computing the approximated solution of Riccati integral equations. In [2], the Riccati solution is applied in a sensor placement problem, which computes optimal sensor locations that minimize the trace of the covariance operator of the Kalman filter of a diffusion-advection process. The same cost has been used in an optimization framework for mobile sensor’s motion planning in [1]. The existence of a solution of the optimization problem is established under the assumption that the integral kernel of the output operator is continuous with respect to the location of the sensor [1, Definition 4.5]. This assumption permits the construction of a compact set of operators [1, Lemma 4.6] over which the cost function is continuous, and hence establishes the existence of a solution. The assumption is also made on the input operator in this paper, which allows the derivation of a vital result on the Riccati operator’s continuity with respect to the actuator trajectory (Lemma 3). The continuity property plays a crucial role in establishing the existence of the proposed problem’s solution and the convergence to the exact optimal solution of the approximate optimal solution. The existence of a solution is established in using the fact that a weakly sequentially lower semicontinuous function attains its minimum on a weakly sequentially compact set over a normed linear space. In addition to the assumptions made for the existence of a solution, a stringent (yet with reasonable physical interpretation) requirement is placed on the admissible set to yield compactness, which leads to convergence of the approximate optimal solution. The convergence is in the sense that when evaluating the exact and approximate optimal solutions by the original cost function, the difference becomes arbitrarily small as the dimension of approximation increases.

The contributions of this paper are threefold. First, we propose an optimization framework for controlling a PDE-modeled system using a team of mobile actuators. The framework incorporates both controlling the process and steering the mobile actuators. The former is handled by the linear-quadratic regulator of the PDE. The latter is taken care of by designing generic cost functions that address the constraints and limitations of the vehicles carrying the actuators. Second, existence conditions of a solution of the proposed problem are established. It turns out that the conditions are generally satisfied in engineering problems, which allows the results to be applied to a wide range of applications. Third, conditions are also established under which the optimal solution computed using approximations converges to the exact optimal solution. The convergence is in the sense that the cost function of the exact problem evaluated at these two solutions becomes arbitrarily close as the dimensional of the approximation goes to infinity. The convergence is verified in numerical studies and confirms the appropriateness of the optimal solution of the approximation.

The proposed framework is well-suited for the limited onboard resources of mobile actuators in the following two aspects: (1) it adopts a finite-horizon optimization scheme that characterizes the resource limitation more precisely than the alternative approaches that do not specify a terminal time, such as an infinite-horizon optimization or Lyapunov-based method; and (2) it provides an intermediate step for the optimization problem that characterizes the limited resources as inequality constraints, because the constraints can be used to augment the cost function and turned into the proposed form using the method of Lagrange multipliers. Potential applications of this work include forest firefighting using unmanned aerial vehicles and oil spill removal or harmful algae containment using autonomous skimmer boats. A preliminary version of this paper [4] considered controlling a 1D diffusion process by a team of mobile actuators. The results in this paper extends the controlled process in [4] to a 2D diffusion-advection process and generalizes the mobility cost therein. Furthermore, the proofs of the existence of a solution and convergence of the approximated optimal solution are presented for the first time in this paper. The results for a dual estimation framework can be found in [5].

The paper adopts the following notation. The symbols $\mathbb{R}$, $\mathbb{R}^{+}$, and $\mathbb{N}$ denote the set of real numbers, nonnegative real numbers, and nonnegative integers, respectively. The boundary of a set $M$ is denoted by $\partial M$. The $n$-nary Cartesian power of a set $M$ is denoted by $M^{n}$. A continuous embedding is denoted by $\hookrightarrow$. We use $\left|\cdot\right|$ and $\left\lVert\cdot\right\rVert$ for the norm defined on a finite- and infinite-dimensional space, respectively, with subscript indicating type. The superscript ^∗ denotes an optimal variable or an optimal value, whereas ^⋆ denotes the adjoint of a linear operator. The transpose of a matrix $A$ is denoted by $A^{\top}$. An $n\times n$-dimensional identity matrix is denoted by $I_{n}$. We denote by $0_{n\times m}$ and $1_{n\times m}$ an $n\times m$-dimensional matrix with all entries being $0$ and $1$, respectively. The term guidance refers to the steering of the mobile actuators, whereas the term control refers to the control input to the DPS. For an optimization problem (P0) that minimizes cost function $J(\cdot)$ over variable $x$ subject to constraints, we use $J_{\text{(P0)}}(x)$ to denote the cost function of (P0) evaluated at $x$. Specifically, $J^{*}_{\text{(P0)}}(x)$ indicates that the optimal value of (P0) is attained when the cost function is evaluated at $x$.

Section II introduces relevant mathematical background, including representation of a diffusion-advection equation by an infinite-dimensional system, the associated LQ optimal control, and its finite-dimensional approximation. Section III introduces the proposed optimization problem and its equivalent problem. Conditions for the existence of a solution are stated. Section IV details the computation of an optimal solution using finite-dimensional approximations. Conditions for the convergence to the exact optimal solution of the approximate optimal solution are stated. A gradient-based method is applied to find an optimal solution. Section V provides two numerical examples to illustrate optimal guidance and control solved by the proposed method. Section VI summarizes the paper and discusses ongoing work.

This paper is motivated by the problem of controlling the following diffusion-advection process on a two-dimensional spatial domain $\Omega=[0,1]\times[0,1]$ with a team of $m_{a}$ mobile actuators:

where $z(\cdot,\cdot,t)$ is the state at time $t$, $\mathbf{v}\in\mathbb{R}^{2}$ is the velocity field for advection, and $u_{i}$ is the control implemented by actuator $i$, with the actuation characterized spatially by $\mathcal{B}_{i}$. The state $z$ lives in the state space $L^{2}(\Omega)$. A representative model of the actuation dispensed by each actuator is Gaussian-shaped and centered at the actuator $i$’s location $(x_{i},y_{i})$ with a bounded support such that

where the parameter $\sigma_{i}$ determines the spatial influence of the actuation, which is concentrated mostly at the location of the actuator and disperses to the surrounding with an exponential decay.

This paper seeks to derive the guidance and control input of each actuator such that the state $\mathcal{Z}$ of the abstract linear system (7) can be driven to zero. Specifically, consider the following problem:

where $J_{\text{m}}(\xi,p):=\int_{0}^{t_{f}}h(\xi(t),t)+g(p(t),t)\text{d}t+h_{f}(\xi(t_{f}))$ is the cost associated with the motion of the actuators, named the mobility cost, such that the mappings $h:\mathbb{R}^{n}\times[0,t_{f}]\rightarrow\mathbb{R}^{+}$ and $g:\mathbb{R}^{m}\times[0,t_{f}]\rightarrow\mathbb{R}^{+}$ evaluate the running state cost and running guidance cost, respectively, and the mapping $h_{f}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{+}$ evaluates the terminal state cost.

The running state cost $h(\cdot,\cdot)$ may characterize restrictions to actuator state. For example, a Gaussian-type function with its peak in the center of the spatial domain, i.e.,

where $\sigma_{x}(t),\sigma_{y}(t)>0$ and $x,y\in[0,1]$, can model a hazardous field that may shorten the life span of an actuator. The integral of this function in the interval $[0,t_{f}]$ evaluates the accumulated exposure of the mobile actuator along its trajectory, which may need to be contained as small as possible (see [11]). Another example is the artificial potential field [16], cast as a soft constraint, that penalizes the trajectory when it passes an inaccessible region such as an obstacle. The running guidance cost $g(\cdot,\cdot)$ may be the absolute value or a quadratic function of the guidance, which characterizes the total amount (of fuel) or energy for steering, respectively. And the terminal state cost $h_{f}(\cdot)$ may characterize restrictions of the terminal state of the mobile actuators. For example, if an application specifies terminal positions, then $h_{f}(\cdot)$ may be a quadratic function that penalizes the deviation of the actual terminal positions.

The formulation in (P) provides an intermediate step for minimizing the PDE cost subject to mobility constraints, in addition to the dynamics constraints. The mobility constraints are characterized by inequalities of $h_{f}(\cdot)$ and the integrals of $h(\cdot,\cdot)$ and $g(\cdot,\cdot)$, because these constraints can be used to augment the cost function and turned into the form of (P) using the method of Lagrange multipliers.

An equivalent problem of (P) can be derived using Lemma 2. For an arbitrary admissible guidance $p$, the actuator trajectory $\xi$ is determined following the dynamics (6), which also determines the input operator $\mathcal{B}(\xi(\cdot),\cdot)$. By Lemma 2, the control $u$ that minimizes the cost function of (P)—specifically, the PDE cost $J(\mathcal{Z},u)$—is given by (10), and the minimum PDE cost is $\langle\mathcal{Z}_{0},\Pi(0)\mathcal{Z}_{0}\rangle$, where $\Pi(0)$ is the mild solution of (12) with actuator trajectory steered by guidance $p$. Hence, we derive the following problem equivalent to (P):

where $\Pi(0)$ is defined in (12) with $t=0$.

To prove the existence of a solution to (P1), we make the following assumptions on the admissible set of guidance and the functions composing the mobility cost:

1. (A9)

   The set of admissible guidance $P\subset\mathbb{R}^{m}$ is closed and convex.

2. (A10)

   The mappings $h:\mathbb{R}^{n}\times[0,t_{f}]\rightarrow\mathbb{R}^{+}$, $g:\mathbb{R}^{m}\times[0,t_{f}]\rightarrow\mathbb{R}^{+}$, and $h_{f}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{+}$ are continuous. For every $t\in[0,t_{f}]$, the function $g(\cdot,t)$ is convex.

3. (A11)

   There exists a constant $d_{1}>0$ with $g(p,t)\geq d_{1}|p|_{2}^{2}$ for all $(p,t)\in P\times[0,t_{f}]$.

Assumptions (A9)–(A11) are generally satisfied in applications with vehicles carrying the actuators. Assumption (A9) is a physically reasonable characterization of the steering of a vehicle, where the admissible steering is generally a continuum with attainable limits within its range. Assumption (A10) places a general continuity requirement on the cost functions and a convexity requirement on the steering cost function. Assumption (A11) requires the function $g(p,t)$ to be bounded below by a quadratic function of the guidance $p$ for all $t$, which is generally satisfied, e.g., with $g$ itself being a quadratic function of $p$. These assumptions are applied in Theorem 6 below regarding the existence of a solution of (P1), whose proof is in Appendix A. Subsequently, the solution to (P1) can be used to reconstruct the solutions to (P), which is stated in Theorem 7 with its proof in Appendix B.

The equivalent problem (P1) allows us to search for an optimal guidance $p$ such that the mobility cost plus the optimal PDE cost is minimized. The control is no longer an optimization variable, because it is determined by the LQR of the abstract linear system for arbitrary trajectories of the mobile actuators.

Approximation of the infinite-dimensional terms in problem (P) is necessary when computing the optimal control and guidance. Hence, we replace the PDE cost and dynamics of (P) by (16) and (13), respectively, and obtain the following approximate problem (AP):

Similar to (P), problem (AP) can be turned into an equivalent form using LQR results for a finite-dimensional system:

where $\Pi_{N}(0)$ is defined in (15) with $t=0$. Analogous to Theorems 6 and 7, the existence of a solution of (AP1) and how to use its solution to reconstruct a solution for (AP) are stated in Theorem 8 below, whose proof is presented in Appendix C.

An extension to Theorem 8 is that an optimal feedback control can be obtained from (17) whenever the optimal guidance is solved from (AP) or (AP1). Basically, when the trajectory is determined via the optimal guidance, a feedback control can be implemented.

To establish convergence to the solution of (P1) of (AP1)’s solution, we need to restrict the set of admissible guidance to a smaller set as introduced below in assumption (A12).

1. (A12)

   There exist $p_{\max}>0$ and $a_{\max}>0$ such that the set of admissible guidance is $\mathcal{P}(p_{\max},a_{\max}):=\{p\in C([0,t_{f}];P):|p(t)|$ is uniformly bounded by $p_{\max}$ and $|p(t_{1})-p(t_{2})|\leq a_{\max}|t_{1}-t_{2}|,\ \forall t_{1},t_{2}\in[0,t_{f}]\}$.

There are two perspectives to interpreting the assumption (A12). Mathematically, (A12) requires the admissible guidance to be a continuous function that is uniformly bounded and uniformly equicontinuous. These two properties yield the sequential compactness of the set $\mathcal{P}(p_{\max},a_{\max})$ by the Arzelà-Ascoli Theorem [30]. Practically, (A12) requires the input signal to be continuous and have bounds $p_{\max}$ and $a_{\max}$ on the magnitude and the rate of change, respectively. This requirement is reasonable and checkable because a continuous signal is commonly used for smooth operation, and the bounds on magnitude and changing rate are due to the physical limits of the motion of a platform. For example, in the case of single integrator dynamics where $p$ is the velocity command, $p_{\max}$ and $a_{\max}$ refer to the maximum speed and maximum acceleration, respectively. Moreover, since time discretization of the signal is applied when computing the optimal guidance, as long as the bound $p_{\max}$ on the magnitude of the signal is determined, then the changing rate is bounded by $a_{\max}:=2p_{\max}/\Delta t_{\min}$ for the smallest discrete interval length $\Delta t_{\min}$. Theorem 9 below states the convergence of the approximate optimal solution with its proof in Appendix D.

The convergence stated in Theorem 9 is established based on several earlier stated results, including

1. 1.

   the input operator’s continuity with respect to location (assumption (A4)), which leads to the continuity of the PDE cost with respect to actuator trajectory (Lemma 3);

2. 2.

   existence of the Riccati operator (Lemma 1) and convergence of its approximation (Theorem 4); and

3. 3.

   sequential compactness of the set of admissible guidance (assumption (A12)), which leads to the continuity of the cost function with respect to guidance (Lemma 12).

Notice that these key results, in an analogous manner, are also required in [24] when establishing the convergence to the exact optimal actuator locations of the approximate optimal locations [24, Theorem 3.5], i.e.,

1. 1.

   continuity with respect to location and compactness of the input operator [24, Theorem 2.6], which lead to continuity of the Riccati operator with respect to actuator locations [24, Theorem 2.6];

2. 2.

   existence of the Riccati operator [24, Theorem 2.3] and the convergence of its approximation [24, Theorem 3.1]; and

3. 3.

   sequential compactness of the set of admissible locations, which is inherited from the setting that the spatial domain is closed and bounded in a finite-dimensional space.

Although the establishment of convergence is similar to the one in [24], the cost function and type of Riccati equation are different: we have the quadratic PDE cost plus generic mobility cost and differential Riccati equation in this paper for control and actuator guidance versus the Riccati operator’s norm as cost function and algebraic Riccati equation in [24] for actuator placement. The similarity comes from the infinite-dimensional nature of PDEs such that approximation is necessary for computation, and convergence in approximation qualifies the approximate optimal solutions.