Finite-Horizon Optimal Control of Boolean Control Networks: A Unified Graph-Theoretical Approach

TL;DR

This paper presents a graph-theoretical framework for finite-horizon optimal control of Boolean control networks, enabling efficient solutions for complex problems including time-variant costs and various constraints, validated on a genetic network example.

Contribution

It introduces a unified graph-based approach for FHOC of BCNs, capable of handling general problems with time-variant costs and improving computational efficiency over existing methods.

Findings

Developed weighted state transition graphs and algorithms for shortest path problems.

Proved equivalence between FHOC and shortest path problems in specific graphs.

Validated approach on E. coli genetic network, showing significant time reduction.

Abstract

This paper investigates the finite-horizon optimal control (FHOC) problem of Boolean control networks (BCNs) from a graph theory perspective. We first formulate two general problems to unify various special cases studied in the literature: (i) the horizon length is $\textit{a priori}$ fixed; (ii) the horizon length is unspecified but finite for given destination states. Notably, both problems can incorporate time-variant costs, which are rarely considered in existing work, and a variety of constraints. The existence of an optimal control sequence is analyzed under mild assumptions. Motivated by BCNs' finite state space and control space, we approach the two general problems in an intuitive and efficient way under a graph-theoretical framework. A weighted state transition graph and its time-expanded variants are developed, and the equivalence between the FHOC problem and the shortest path problem in specific graphs is established rigorously. Two custom algorithms are developed to find the shortest path and construct the optimal control sequence with lower time complexity, though technically a classical shortest-path algorithm in graph theory is sufficient for all problems. Compared with existing algebraic methods, our graph-theoretical approach can achieve state-of-the-art time efficiency while targeting the most general problems. Furthermore, our approach is the first one capable of solving Problem (ii) with time-variant costs. Finally, the Ara operon genetic network in $\textit{E. coli}$ is used as a benchmark example to validate the effectiveness of our approach, and the results of two tasks show that our approach can dramatically reduce the running time.

Tables

Table 1. TABLE I: Time complexity comparison between existing work and the proposed approach on finite-horizon optimal control

Problem type	Task characteristics	Time complexity
		Existing work	Proposed approach
Problem 1	Mayer-type optimal control (only terminal cost)	$O\left(TM{N}^2\right)$ [12, 13]	$O\left(TMN\right)$ (Algorithm 2)
	Minimum-energy control (only time-invariant stage cost)	$O\left(T{N}^3\right)$ [14, Algorithm 3.2]
	Time-discounted stage cost (no terminal cost)	$O\left(MN+N^3{\log }_{2}T\right)$ [17, Theorem 3] $O\left(M^2{N}^3+TMN\left(N+M\right)\right)$ [18, Proposition 4.1]
	Both terminal cost and stage cost (general form, time-invariant)	$O\left(MN+N^3{\log }_{2}T\right)$ [17, Corollary 2] $O\left(TMN+TN\log \left(TN\right)\right)$   [23, Theorem 2.14]
	Both terminal cost and stage cost (general form, possibly time-variant)	$O\left(TMN\right)$ [16]
Problem 2	Time-optimal control (only constant stage cost equal to 1)	$O\left(M{N}^3\right)$$†$   [15, 32]	$O\left(MN+N\log N\right)$ (Algorithm 3)
	Minimum-energy control (time-invariant stage cost, time not fixed)	$O\left(N^4\right)$ [14, Algorithm 3.3 ]
	Time-invariant stage cost (general form, no terminal cost)	$O\left(MN+N\log N\right)$   [23, Theorem 2.7]
	Time-variant stage cost and (or) terminal cost	—	$O\left(N^2\left(M+2\log N\right)\right)$ (Algorithm LABEL:S-alg:_dijkstra_case2)

Table 2. TABLE II: Boolean Model of the Ara Operon Network

Target node	Boolean update rule
$A$	$A_e\wedge T$
$A_m$	$\left(A_{em}\wedge T\right)\vee A_e$
$A_{r{a}_{+}}$	$\left(A_m\vee A\right)\wedge A_{r{a}_{-}}$
$C$	$\neg G_e$
$E$	$M_s$
$D$	$\neg A_{r{a}_{+}}\wedge A_{r{a}_{-}}$
$M_S$	$A_{r{a}_{+}}\wedge C\wedge \neg D$
$M_T$	$A_{r{a}_{+}}\wedge C$
$T$	$M_T$

Table 3. TABLE III: Running Time of Various Methods for Optimal Control of the Ara Operon Network

Problem	Task 1					Task 2
Method	[14, Algorithm 3.2]	[16]	[17]	[23, Theorem 2.14]	Ours	[14, Algorithm 3.3]	[15]	[23, Theorem 2.7]	Ours
Running time (s)	350.43	0.23	402.79	0.17	0.10	19651.33	636.72	0.01	0.01

Equations

$$A⋉B=(A\otimes I_{\frac{s}{n}})(B\otimes I_{\frac{s}{p}}),$$

$$f(x_1,x_2,\cdots ,x_n)=M_f{x}_1{x}_2\cdots x_n.$$

$$\Sigma :\{\begin{array}{l}{x}_1(t+1)=f_1(x_1(t),\cdots ,x_n(t),u_1(t),\cdots ,u_m(t))\\ x_2(t+1)=f_2(x_1(t),\cdots ,x_n(t),u_1(t),\cdots ,u_m(t))\\ ⋮\\ x_n(t+1)=f_n(x_1(t),\cdots ,x_n(t),u_1(t),\cdots ,u_m(t)),\end{array}$$

$$x(t+1)=Lu(t)x(t),$$

$$w(p)=\sum _{j=0}^{k-1}w(v_{i_j},v_{i_{j+1}}).$$

$$\underset{u}{\min }{J}_{\text{T}}(u)=h_{\text{T}}(x(T))+\sum _{t=0}^{T-1}g(x(t),u(t),t),$$

$$\text{s.t.}\{\begin{array}{l}x(t+1)=Lu(t)x(t)\\ x(t)\in C_x\\ u(t)\in C_u(x(t))\\ x(0)=x_0\\ x(T)\in \Omega \end{array},$$

$$h_{\text{T}}(x(T))=x(T{)}^{\top }{Q}_{h}x(T),$$

$$g(x(t),u(t),t)=\left[\begin{array}{c}x(t{)}^{\top }u(t{)}^{\top }\end{array}\right]\left[\begin{array}{cc}Q& S\\ S^{\top }& R\end{array}\right]\left[\begin{array}{c}x(t)\\ u(t)\end{array}\right],$$

$$\underset{u}{\min }J(u)=h(x(K),K)+\sum _{t=0}^{K-1}g(x(t),u(t),t),$$

$$\text{s.t.}\{\begin{array}{l}x(t+1)=Lu(t)x(t)\\ x(t)\in C_x\\ u(t)\in C_u(x(t))\\ x(0)=x_0\\ x(K)\in \Omega \end{array},$$

$${\mathcal{R}}_1({\delta }_N^i)=\{{\text{Col}}_i({\text{Blk}}_k(L))|k\in [1,M]\},$$

$${\mathcal{R}}_1({\delta }_N^i)=\{{\text{Col}}_i({\text{Blk}}_k(L))|k\in [1,M],{\delta }_M^k\in C_u({\delta }_N^i)\}\cap C_x.$$

$$U^{ij}=\{{\delta }_M^k\in C_u({\delta }_N^i)|k\in [1,M],{\text{Col}}_i({\text{Blk}}_k(L))={\delta }_N^j\}.$$

$${\mathcal{R}}_d({\delta }_N^i)={\mathcal{R}}_1({\mathcal{R}}_{d-1}({\delta }_N^i)),d\ge 2.$$

$$E=\{({\delta }_N^i,{\delta }_N^j)|{\delta }_N^i\in V,{\delta }_N^j\in {\mathcal{R}}_1({\delta }_N^i)\},$$

$${\Sigma }_1:\{\begin{array}{l}{x}_1(t+1)=x_2(t)\wedge (u_1(t)\oplus x_3(t))\\ x_2(t+1)=\neg x_1(t)\\ x_3(t+1)=u_2(t)\oplus x_2(t)\end{array}.$$

$$\{\begin{array}{ll}{C}_x& ={\Delta }_8\setminus \{{\delta }_8^8\}\\ C_u(x)& =\{\begin{array}{l}\{{\delta }_4^3,{\delta }_4^4\},\text{if }x={\delta }_8^6\\ {\Delta }_4\setminus \{{\delta }_4^2\},\text{otherwise}\end{array}\end{array}$$

$$w({\delta }_{N,t}^i,{\delta }_{N,t+1}^j)=\{\begin{array}{l}{\min }_{{\delta }_M^k\in U^{ij}}g({\delta }_N^i,{\delta }_M^k,t),t\in [0,T-1]\\ h_{\text{T}}({\delta }_N^i),t=T\end{array}$$

$$u_t^{ij}=\underset{{\delta }_M^k\in U^{ij}}{\mathrm{argmin}}g({\delta }_N^i,{\delta }_M^k,t),t\in [0,T-1].$$

$$g(x(t),u(t),t)=u(t{)}^{\top }Qu(t)+t,h_{\text{T}}(x(T))=x(T{)}^{\top }Rx(T),$$

$$p=\big{<}\delta_{N,0}^{i_{0}},\delta_{N,1}^{i_{1}},\cdots,\delta_{N,T}^{i_{T}},\delta_{N,T+1}^{0}\big{>},$$

$$w(p)=h_{\text{T}}({\delta }_N^{i_T})+\sum _{t=0}^{T-1}g({\delta }_N^{i_t},u_t^{i_t{i}_{t+1}},t).$$

$$J_{\text{T}}(u)=h_{\text{T}}({\delta }_N^{i_T})+\sum _{t=0}^{T-1}g({\delta }_N^{i_t},{\delta }_M^{k_t},t).$$

$$p^{*}=\big{<}\delta_{N,0}^{i_{0}},\delta_{N,1}^{i_{1}^{*}},\cdots,\delta_{N,T}^{i_{T}^{*}},\delta_{N,T+1}^{0}\big{>},$$

$$u^{\ast }=\{u_0^{i_0{i}_1^{\ast }},u_1^{i_1^{\ast }{i}_2^{\ast }},\cdots ,u_{T-1}^{i_{T-1}^{\ast }{i}_T^{\ast }}\}.$$

$$P({\delta }_{N,t}^j)=\{{\delta }_{N,t-1}^i|({\delta }_{N,t-1}^i,{\delta }_{N,t}^j)\in E\},t\in [1,T+1]$$

$$F({\delta }_{N,t}^j)=\underset{{\delta }_{N,t-1}^i\in P({\delta }_{N,t}^j)}{\min }F({\delta }_{N,t-1}^i)+w({\delta }_{N,t-1}^i,{\delta }_{N,t}^j),$$

$$J(u)$$

$$\displaystyle+\sum_{t=j}^{k-1}\left[g(x(t),u(t),t)-g\big{(}x(t),u(t),t-r\big{)}\right].$$

$$w({\delta }_N^i,{\delta }_N^j)=\underset{{\delta }_M^k\in U^{ij}}{\min }g({\delta }_N^i,{\delta }_M^k,\alpha ),({\delta }_N^i,{\delta }_N^j)\in E.$$

$$V\leftarrow V\cup \{{\delta }_N^0\},E\leftarrow E\cup E^0$$

Full text

I Introduction

Boolean network (BN) was first proposed by Kauffman [1] to model gene regulatory networks, where each gene is assigned a Boolean variable to represent its expression state. The BN model has thereafter attracted increasing research interest in various fields, including studies on biomolecular networks in systems biology [2], therapeutic interventions in clinical treatment [3], and the contagion dynamics during a financial crisis [4], just to name a few. In a BN, the binary variables interact with each other through Boolean functions, and exogenous (binary) inputs can be injected into these functions to affect the network dynamics, which is commonly referred to as a Boolean control network (BCN) [5]. In this study, we focus on finite-horizon optimal control (FHOC) of BCNs, which is useful for optimal therapeutic intervention strategy design in medical applications [3]. Note that we view the finite horizon in a more general sense: either a fixed horizon length, or an unknown but finite horizon length with respect to a given destination state. These two types of FHOC problems are referred to as fixed-time optimal control and fixed-destination optimal control respectively in this paper.

In the last decade, a new matrix product called the semi-tensor product (STP), which can convert a BN (BCN) into an algebraic state-space representation (ASSR), has been developed by Daizhan Cheng et al. [6, 7]. The ASSR provides a systematic framework to study a wide range of control-theoretical problems related to BCNs, such as their controllability [5], observability [5, 8], stabilization [9], and various controller synthesis problems [10, 11], among others. A number of optimal control problems for BCNs have been investigated in recent years using the STP and ASSR tools as well. In [12], the Mayer-type optimal control problem (i.e., only terminal cost is considered) for single-input BCNs is addressed, and a necessary condition analogous to the Pontryagin’s maximum principle is derived, which has been later extended to multi-input BCNs [13]. The minimum-energy control and minimum-time control of BCNs are investigated in [14] and [15] respectively. More general FHOC problems involving both stage cost and terminal cost are considered in [16], and the solution is given by a recursive algorithm as an analogy to the difference Riccati equation for discrete-time linear systems. More recently, Ref. [17] targets the time-discounted stage cost and introduces a recursive algorithm based on a data structure called the optimal input-state transfer graph. The same problem is also investigated in [18], and a recursive solution for receding horizon optimal control of mix-valued probabilistic logical networks is obtained. In parallel to the study of FHOC problems, the more challenging infinite-horizon counterparts have also been attempted recently in several contributions using STP-based algebraic methods. For example, infinite-horizon optimal control with average cost is studied in [19, 20, 16], and the time-discounted cost case is examined in [21, 17, 22].

While we appreciate the above successes achieved with algebraic methods powered by the STP theory, one issue is that distinct methods are developed to solve different problems as reviewed above, even if we only consider FHOC. It appears that most existing work only deals with certain special cases of FHOC, for example, the minimum-energy control [14], the minimum-time control [15], the Mayer-type problem [12], the two kinds of Lagrange-type problems [23], and the time-discounted finite-horizon problem [18], though they may share a lot in common. As mentioned in [16] and [24], there are many other types of FHOC problems as well as lots of practical limitations, for instance, optimal control subject to various constraints [25, 24], and more challenging problems with general time-variant costs. Consequently, a versatile approach instead of fragmented methods to solve all these common problems is highly desirable, which partly motivates this research. One goal of our study is to unify various tasks into an integrated framework and to develop systematic algorithms to handle the most general problems.

One serious concern about the algebraic approaches in BCN studies is the high computational complexity, since the number of operations grows exponentially with respect to the network size [5, 7, 11]. As a result, these approaches may quickly become computationally intractable as the number of state variables increase. In fact, most control-theoretical problems of BCNs are NP-hard [15] due to the combinatorial nature of Boolean state variables, i.e., a BCN with $n$ variables has up to $N\vcentcolon=2^{n}$ states, a phenomenon well-known as the curse of dimensionality, such as NP-hardness of controllability [26] and observability [8]. However, this doesn’t mean we are hopeless: many algebraic approaches run in a high-order polynomial time in terms of $N$ (see Table I for more details), which still leaves vast room for improvement by designing more efficient algorithms to decrease the order of the polynomial [27, 9]. This forms another objective of this research: we aim to reduce the computational complexity of FHOC for BCNs.

To pursue more efficient algorithms, we notice that a BCN is characterized by its finite state space, finite control space, and deterministic state transitions. A direct consequence of this property is that the dynamics of a BCN can be adequately described by a state transition graph (STG). Going further, methods originating from graph theory appear to be promising for investigations of BCNs. For example, computationally efficient methods for controllability [27, 28] and stabilization [9] of BCNs all resort to certain graph-theoretical algorithms. Regarding the optimal control of BCNs, a couple of pioneering studies exist that attempt to combine the ASSR with tools from graph theory to improve computational efficiency. A typical example is the employment of Floyd-like algorithms, inspired by the Floyd-Warshall algorithm in graph theory, in both finite-horizon minimum-energy control [14] and infinite-horizon problems [20, 21]. An update-to-date study [23] handles two kinds of Lagrange-type optimal control problems using the Dijkstra’s algorithm instead. The above two algorithms are both initially developed to find shortest paths in a weighted graph. One limitation of both [14] and [23] is that they only consider time-invariant costs and a special class of FHOC problems. Motivated by these pioneering work borrowing tools from graph theory, we attempt to advance further and unify all common FHOC problems into an elegant and efficient graph-theoretical framework.

The contributions of this paper are listed in three folds. First, we unify all common types of FHOC problems, including time-variant costs and various constraints, into two general problems depending on whether the horizon length is prespecified, which are subsequently reduced to shortest path problems by constructing particular variants of the STG. The correctness of this reduction is proved rigorously. Existence of an optimal control sequence is analyzed under mild conditions for both problems. Second, we develop two intuitive algorithms to solve the above shortest path problems with superior efficiency. To be specific, only one algebraic method can achieve the same worst-case time complexity as ours, but our approach still tends to have better practical performance, which is demonstrated by two optimal control tasks for a genetic network in the bacteria E. coli. Third, as far as we know, there are currently no published results on the fixed-destination optimal control problem with time-variant costs, and our approach can handle it effectively using the identical graph-theoretical methodology.

The remainder of this paper is organized as follows. In Section II, we present some background knowledge about the STP, the ASSR, and the shortest path problem in graph theory. Section III introduces two general problems that can incorporate all specific FHOC problems studied in the current literature. In Section IV, the construction of the STG from a given initial state is discussed. After that, we detail the equivalence between the two FHOC problems and the shortest path problem in dedicated graphs and propose two efficient algorithms to solve them in Section V and VI respectively. The time complexity of our approach is compared with that of existent work in Section VII, and the practical running time of various methods to complete two tasks for the Ara operon network of E. coli is measured in Section VIII. Finally, Section IX gives some concluding remarks.

II Preliminaries

II-A Notations

For statement ease, the following notations [7, 6] are used.

$|\mathcal{S}|$ denotes the size (i.e., cardinality) of a set $\mathcal{S}$. 2. 2.

Let $\mathbb{R}$ and $\mathbb{N}$ denote the set of real numbers and nonnegative integers respectively. $[l,r]\vcentcolon=\{l,l+1,\cdots,r-1,r\}$. 3. 3.

$f\geq B$ means a function $f$ is bounded below by $B$. 4. 4.

$\mathcal{M}_{p\times q}$ denotes the set of all $p\times q$ matrices. Given $A\in$ 5. 5.

$\textrm{Col}_{i}(M)$ denotes $i$-th column of a matrix $M$, and $M_{ij}$ denotes the $(i,j)$-th entry of $M$. 6. 6.

$\delta_{n}^{i}\vcentcolon=\textrm{Col}_{i}(I_{n})$, where $I_{n}\in\mathcal{M}_{n\times n}$ is the identity matrix. $\Delta_{n}\vcentcolon=\{\delta_{n}^{i}|i=1,2,\cdots,n\}$, and $\Delta\vcentcolon=\Delta_{2}$. 7. 7.

A matrix $M=[\delta_{n}^{i_{1}}\;\delta_{n}^{i_{2}}\;\cdots\;\delta_{n}^{i_{q}}]\in\mathcal{M}_{n\times q}$ with $\delta_{n}^{i_{k}}\in\Delta_{n},\forall k\in[1,q]$, is called a logical matrix. Let $\mathcal{L}_{n\times q}$ denote the set of all $n\times q$ logical matrices. 8. 8.

A matrix $A\in\mathcal{M}_{n\times mn}$ can be rewritten into a block form $A=[\text{Blk}_{1}(A)\;\text{Blk}_{2}(A)\;\cdots\;\text{Blk}_{m}(A)]$, where $\text{Blk}_{i}(A)\in\mathcal{M}_{n\times n}$ is the $i$-th square block of $A$. 9. 9.

Logical operators [6]: $\land$ for conjunction, $\lor$ for disjunction, $\lnot$ for negation, and $\oplus$ for exclusive disjunction.

II-B STP of Matrices and ASSR of BCNs

This section revisits some necessary background knowledge about the STP and the ASSR of BCNs developed by Daizhan Cheng et al. [6, 7, 5].

Definition 1

[19]** The STP of two matrices $A\in\mathcal{M}_{m\times n}$ and $B\in\mathcal{M}_{p\times q}$ is defined by

[TABLE]

where $\otimes$ denotes the Kronecker product, and $s$ is the least common multiple of $n$ and $p$.

Remark 1

The STP is essentially a generalization of the standard matrix product, and all major properties of the standard matrix product remain valid under STP [19]. Thus, we will omit the symbol $\ltimes$ in the remainder when no confusion is caused. That is, all matrix products refer to STP by default.

A logical function can be conveniently expressed in a multilinear form via STP. In this form, a Boolean value is identified by a vector as $\text{TRUE}\sim\delta_{2}^{1}$ and $\text{FALSE}\sim\delta_{2}^{2}$.

Lemma 1

[6]** Given any Boolean function $f(x_{1},x_{2},\cdots,x_{n}):\Delta^{n}\rightarrow\Delta$, there exists a unique matrix $M_{f}\in\mathcal{L}_{2\times 2^{n}}$, called the structure matrix, such that

[TABLE]

Interested readers may refer to [6] for computation of $M_{f}$. For simplicity, set $\ltimes_{i=1}^{n}A_{i}\vcentcolon=A_{1}\ltimes\cdots\ltimes A_{n}$.

A general BCN $\Sigma$ with $n$ state variables (i.e., $n$ nodes in the network) and $m$ control inputs can be described as

[TABLE]

where $x_{i}(t)\in\Delta$ denotes the value of the $i$-th variable at time $t$, and $f_{i}:\Delta^{n+m}\rightarrow\Delta$ is the $i$-th Boolean function, $i\in[1,n]$. Additionally, $u_{j}(t)\in\Delta$ denotes the $j$-th control input at time $t$, $j\in[1,m]$. Set $x(t)\vcentcolon=\ltimes_{i=1}^{n}x_{i}(t)$ and $u(t)\vcentcolon=\ltimes_{j=1}^{m}u_{j}(t)$. Clearly, we have $x(t)\in\Delta_{N}$ and $u(t)\in\Delta_{M}$, where $N\vcentcolon=2^{n}$ and $M\vcentcolon=2^{m}$. $N$ and $M$ will be used through the whole text.

The ASSR of the BCN in (2) is given by [6, 14]:

[TABLE]

where $L\in\mathcal{L}_{N\times MN}$, named the network transition matrix, is computed by $\text{Col}_{j}(L)=\ltimes_{i=1}^{n}\text{Col}_{j}(M_{f_{i}}),\forall j\in[1,MN]$, where $M_{f_{i}}\in\mathcal{L}_{2\times MN}$ is the structure matrix of $f_{i}$ in (2).

II-C Shortest Path Problem

The core of our graph-theoretical approach for FHOC is to transform the original problem into a shortest path problem in a particular graph and then locate the shortest path efficiently. We brief the shortest path problem in graph theory below.

Given a directed graph $G=(V,E)$, where $V=\{v_{1},v_{2},\cdots,v_{n}\}$ is a set of vertices, and $E=\{(v_{i},v_{j})|v_{i},v_{j}\in V\}$ is a set of directed edges, we assign each edge a real value, called its weight or cost. Denote the weight of the edge from $v_{i}$ to $v_{j}$ by $w(v_{i},v_{j})$. Such a graph is known as a weighted directed graph. A path from a source vertex $v_{i_{0}}\in V$ to a destination vertex $v_{i_{k}}\in V$ is a sequence of vertices connected by edges, denoted by $p=\left<v_{i_{0}},v_{i_{1}},\cdots,v_{i_{k}}\right>$. Let $\epsilon(p)$ and $|p|$ denote the number of edges and vertices in $p$ respectively. If $v_{i_{0}}=v_{i_{k}}$, $p$ is a cycle. $\epsilon(p)$ is also called the length of $p$.

Definition 2

The weight $w(p)$ of a path $p$ is the sum of the weights of its constituent edges:

[TABLE]

A shortest path (SP) from $v_{i_{0}}$ to $v_{i_{k}}$ is any path from $v_{i_{0}}$ to $v_{i_{k}}$ with the minimum weight among all possible paths.

III Problem Formulation

Despite the variety of FHOC problems of BCNs studied in the literature [12, 13, 14, 15, 16, 18, 17, 23], they can essentially be classified into two general types according to whether the horizon length is *a priori *fixed or not. This section will detail the mathematical formulation of both general problems.

III-A Fixed-Time Optimal Control

Roughly speaking, our objective is to construct a fixed-length control sequence for the BCN (3) to optimize a given performance index [25, 16]. The scenario can become much more complicated in practice because of various constraints when designing control strategies. For example, we have to avoid dangerous states in therapeutic intervention [3, 25]. Besides, not all theoretical control inputs are practically realizable, e.g., we may still lack effective means to manipulate specific genes in a GRN, or a medical treatment is possibly unaffordable. Consequently, certain control inputs can be unavailable or just prohibited in specific states [24]. Finally, it is common that we may want to steer the network to a particular state, e.g., to lead a gene regulatory network from a cancerous state (initial state) to a healthy state (terminal state) [29, 30]. Thus, we further enrich the fixed-time optimal control problem by constraining the terminal states.

A general problem reflecting the above intention subject to various constraints is formulated mathematically as follows.

Problem 1

Consider the BCN (3), an initial state $x_{0}\in\Delta_{N}$, and a fixed horizon length $T\in\mathbb{N}$. Fixed-time optimal control of a BCN is to determine an optimal control sequence of length $T$ to the following optimization problem:

[TABLE]

where $u\vcentcolon=(u(0),u(1),\cdots,u(T-1))\in\Delta_{M}^{T}$ is a control sequence; $g:\Delta_{N}\times\Delta_{M}\times\mathbb{N}\rightarrow\mathbb{R}$ is the stage cost function; and $h_{\textrm{T}}:\Delta_{N}\rightarrow\mathbb{R}$ is the terminal cost function. $C_{x}$ denotes the allowed states; $C_{u}$ represents the state-dependent constraints on control inputs; and the terminal set $\Omega\subseteq\Delta_{N}$ denotes the set of desirable terminal states at time $T$.

Remark 2

The horizon length $T$ denotes a finite treatment window in therapeutic applications [3, 25]. Thus, optimal therapeutic intervention strategies can be developed within a treatment window by setting up a proper optimality criterion $J_{\textrm{T}}$ in (1) by domain expert knowledge.

Remark 3

Most existing studies only deal with time-invariant stage cost, that is, the function $g$ in (1) doesn’t really depend on time $t$ (see [20, 15, 23] for examples). Few studies consider a time-dependent $g$ but only in restricted forms such as the time-discounted cost in [17] and [18]. We intend to investigate the most general form (1) directly, where the stage cost $g$ can incorporate time $t$ in any form.

Problem 1 represents a very general setting of fixed-time optimal control problems. A variety of specific finite-horizon problems investigated in existing work can be rewritten easily into this general form. We demonstrate in the sequel how to specialize Problem 1 to some specific fixed-time optimal control problems with different characteristics in the literature.

No terminal state constraints [16, 17, 23]. If the terminal state is not specified, we only need $\Omega=\Delta_{N}$. 2. 2.

A single desired terminal state [15, 14]. Set $\Omega=\{x_{d}\}$, where $x_{d}$ is the single desired terminal state. 3. 3.

Only stage cost [14, 18, 23]. To remove the terminal cost, simply nullify it by $h_{\textrm{T}}(\cdot)\equiv 0$. 4. 4.

Only terminal cost (also known as the Mayer problem) [12, 13]. Just set zero stage costs: $g(\cdot,\cdot,\cdot)\equiv 0$. 5. 5.

Special terminal or stage cost functions [16, 14, 18]. For example, [17] and [18] consider the time-discounted finite-horizon optimal control, which can be expressed by (1) with $g(x(t),u(t),t)=\lambda^{t}c_{g}(u(t),x(t))$, where $0<\lambda<1$ is the discount factor, and $c_{g}(\cdot,\cdot)$ is a time-invariant stage cost. As another example, the energy function in [14] is obtained through $g(x(t),u(t),t)=u^{\top}(t)Qu(t)$, and the more general quadratic cost function in [16] can be implemented straightforwardly by

[TABLE]

where $Q_{h},Q,S,\text{and }R$ are proper weight matrices.

III-B Fixed-Destination Optimal Control

A common task that Problem 1 cannot cover is time-optimal control, which aims to find a control sequence to drive the BCN from a given initial state $x_{0}$ to another destination state $x_{d}$ in minimum time. Such time-optimal control has been widely studied for traditional linear time-invariant (LTI) discrete-time systems, such as the famous deadbeat controller (see [31] for a review). In [15], the minimum-time control of BCNs is first investigated, and the same problem is considered again in [32] but with impulsive disturbances. Such optimality concept can be generalized to other criteria, for example, the minimum-energy control studied in [14], which attempts to steer the BCN to a target state using minimum energy.

Note that although the horizon length is not a priori fixed, it must be finite for a well-posed problem, that is, the specified destination state $x_{d}$ should be reachable from the initial state $x_{0}$ in finite steps. It is the main reason we consider such problems as another class of FHOC problems in a general sense. To further generalize this problem, just like Problem 1, we make it admit a set of destination states and subject to various constraints, which is formalized as follows.

Problem 2

Consider the BCN (3), an initial state $x_{0}\in\Delta_{N}$, and a terminal set $\Omega\subseteq\Delta_{N}$. The fixed-destination optimal control problem of a BCN is to determine an optimal control sequence of a variable length to the optimization problem,

[TABLE]

where $u=(u(0),u(1),\cdots,u(K-1))\in\Delta_{M}^{K}$ represents a control sequence, and $K\in\mathbb{N}$ indicates an unknown but finite horizon length. The terminal cost function $h:\Delta_{N}\times\mathbb{N}\rightarrow\mathbb{R}$ and the stage cost function $g:\Delta_{N}\times\Delta_{M}\times\mathbb{N}\rightarrow\mathbb{R}$ can be time-dependent. $C_{x}$, $C_{u}$, and $\Omega$ denote the state constraint, the control constraint, and desirable destination states respectively.

Remark 4

Despite the apparent similarity in mathematical forms between Problem 1 and Problem 2, they will be treated by distinctly different means to maximize the computational efficiency of each problem. In existing work like [14] and [23], only specific problems with a time-independent stage cost function and a single destination state are investigated, and no terminal cost is considered. However, we argue that it is sensible to set different costs if the desired terminal state is reached at different time, for example, when the desired state refers to a good state of physical health.

Compared with the fixed-time optimal control (Problem 1), there are fewer studies on Problem 2. As far as we know, only the following two custom problems have been investigated in the literature, both with a single destination state, i.e., $\Omega=\{x_{d}\}$, and a time-invariant stage cost function.

Minimum-time control [15, 32]. This kind of control is often referred to as time-optimal control. In Problem 2, simply set $g(\cdot,\cdot,\cdot)\equiv 1$ and $h(\cdot,\cdot)=0$, and we have $J(u)=K$. Thus, we are indeed minimizing the number of steps to steer the BCN (3) from $x_{0}$ to $x_{d}$. 2. 2.

Minimum-energy control [14]. Let $g(x(t),u(t),t)=u^{\top}(t)Qu(t)$ and $h(\cdot,\cdot)=0$, where $Q$ is a positive definite diagonal matrix measuring energy consumption.

IV State Transition Graph

In this section,we introduce the state transition graph (STG) of a BCN, which is useful for both Problem 1 and 2. An efficient algorithm based on bread-first search (BFS) [33] in a graph is developed to construct the STG.

IV-A Reachability of a BCN

To construct the STG, we need to first discuss the reachability of a BCN [5, 7], especially the reachable set of a given intial state, whose definition is given as follows.

Definition 3

The set of states that can be reached from $x(0)=x_{0}\in\Delta_{N}$ at time $t=d$ is $\mathcal{R}_{d}(x_{0})$. Given a set $X\subseteq\Delta_{N}$, let $\mathcal{R}_{d}(X)\vcentcolon=\cup_{x\in X}\mathcal{R}_{d}(x)$ for notational simplicity. The complete reachable set of a given state $x_{0}$ is $\mathcal{R}(x_{0})\vcentcolon=\cup_{d\in\mathbb{N}}\mathcal{R}_{d}(x_{0})$.

$\mathcal{R}(x_{0})$ can be obtained algebraically by calculating powers of the network transition matrix $L$ [5, 7]. However, from a graph-theoretical view, a computationally economical way is to adopt the standard bread-first search (BFS) procedure to build $\mathcal{R}(x_{0})$ iteratively [33, 9]. More interestingly, this method is closely related with the adjacency-list representation of a graph [33], and only requires successive computation of the one-step reachable set $\mathcal{R}_{1}(\cdot)$. Recall that $u(t)\in\Delta_{M}$ and $x(t)\in\Delta_{N}$ are both logical vectors with all zero entries except a single entry of value 1. Thus, given $u(t)=\delta_{M}^{k}$, we have $Lu(t)=\textrm{Blk}_{k}(L)$. The following lemma follows immediately.

Lemma 2

The one-step reachable set of a state $\delta_{N}^{i}$ for a BCN (3) with no constraints is

[TABLE]

and, if there are constraints like Problem 1 (or Problem 2),

[TABLE]

Note that we ignore the terminal constraint in (8) and will handle it later when building specific graphs. It is possible that more than one control input can attain the transition from $\delta_{N}^{i}$ to $\delta_{N}^{j}$. Collect these qualified control inputs into a set $U^{ij}$:

[TABLE]

By Definition 3 and Lemma 2, the $d-$step reachable set of a state $\delta_{N}^{i}$ is obtained efficiently with the recursion:

[TABLE]

IV-B Construction of an STG

As aforementioned, a BCN is characterized by its finite state space, where the finite control inputs coordinate the state transitions deterministically. This feature makes it possible to capture the complete dynamics of a BCN with a directed graph, called the state transition graph (STG).

Definition 4

Consider the BCN (3). Its state transition graph (STG) is a directed graph $G=(V,E)$, where $V=\Delta_{N}$ is the vertex set, i.e., one vertex for each state, and the edge set is

[TABLE]

i.e., one edge for each one-step transition between states.

In this study, we only care about states reachable from an initial state $x_{0}$. We denote such an STG by $G=(V,E,x_{0})$ with $V=\mathcal{R}(x_{0})$. Hereafter we will use the terms vertex and state interchangeably when no ambiguity is caused.

In addition, we may assign a weight to each edge of the STG corresponding to the cost of the each state transition, which will be detailed in following sections.

Remark 5

The input-state graph (in a matrix form), proposed by [21, 7], uses an input-state pair $(\delta_{M}^{i},\delta_{N}^{j})$ as an vertex, leading to $MN$ vertices in total. Our STG with only $N$ or $\mathcal{R}(x_{0})$ vertices is potentially more space and time efficient.

Note from (11) that $\mathcal{R}_{1}(\delta_{N}^{i})$ essentially denotes the successors of the vertex $\delta_{N}^{i}$ in the STG. We thus can build the STG following a BFS procedure and get $\mathcal{R}(x_{0})$ at the same time. BFS is a graph traversal algorithm where the neighbors of a vertex are visited in a FIFO (first-in-first-out) order. Algorithm 1 details the construction of an STG. Running Algorithm 1, we get the reachable set $\mathcal{R}(x_{0})$, i.e., $R$ in the algorithm, and the adjacency list of each vertex, $\mathcal{R}_{1}(x),\forall x\in\mathcal{R}(x_{0})$, which effectively gives the STG $G=(V,E,x_{0})$ by Definition 4.

Time complexity analysis of Algorithm 1. For each state $\delta_{N}^{i}\in\mathcal{R}(x_{0})$, Eq. (7) or (8) needs $M$ operations, and obviously $|\mathcal{R}_{1}(\delta_{N}^{i})|\leq M$. The for loop has at most $M$ iterations; and the while loop runs $|V|$ times because each vertex is enqueued and dequeued exactly once. Thus, the time complexity of Algorithm 1 is $O(M|V|)$, which is equivalent to $O(MN)$ since there are at most $N$ states (vertices). This conforms to the celebrated theorem that BFS runs in linear time with respect to the number of edges and vertices [33], i.e., $O(|V|+|E|)$, because the STG has at most $MN$ edges.

Remark 6

Algorithm 1 implies that, any state $x\in\mathcal{R}(x_{0})$ can be reached from $x_{0}$ in less than $N$ steps, i.e., $\mathcal{R}(x_{0})=\cup_{0\leq d<N}\mathcal{R}_{d}(x_{0})$. This result is intuitive: if a trajectory from $x_{0}$ to $x$ contains more than $N$ states, there must be repetitive ones, and we can remove such cycles to shorten the trajectory.

Example 1

Consider a BCN with $m=2$ control inputs and $n=3$ state variables [14] as follows:

[TABLE]

Its ASSR (3) has a $8\times 32$ transition matrix $L$, which is omitted here to conserve space. For illustration purpose, we set up the state constraints and control constraints arbitrarily below:

[TABLE]

i.e., the state $\delta_{8}^{8}$ is forbidden; the control $\delta_{4}^{2}$ is unavailable to all states; and only control $\{\delta_{4}^{3},\delta_{4}^{4}\}$ is applicable to state $\delta_{8}^{6}$.

The STG $G=(V,E,x_{0})$ of $\Sigma_{1}$ with states reachable from $x_{0}=\delta_{8}^{1}$ subject to the given constraints (13) is obtained by Algorithm 1 and shown in Fig. 1. ∎

V Solve Problem 1 via Dynamic Programming

In this section, we first validate the existence of optimal solutions to Problem 1, and then we introduce a time-expanded variant of the STG and develop an efficient algorithm via dynamic programming to solve Problem 1 based on that graph.

V-A Existence of Optimal Solutions

As we will discuss next, the stage costs and the terminal costs in (1) are mapped to edge weights of a graph to reduce Problem 1 to an SP problem. It is known that the SP problem is well defined only if the graph contains no negative-weight cycles [33]. Nonetheless, since the length of a path in Problem 1 is fixed due to the fixed horizon, this condition is no longer required. Only the following assumption is needed for Problem 1 to avoid an optimal value of negative infinity.

Assumption 1

The stage cost function $g$ and the terminal cost function $h_{\textrm{T}}$ in Problem 1 are both bounded from below.

Now we are ready to present the following conclusion about the existence of optimal solutions to Problem 1.

Proposition 1

Consider Problem 1 under Assumption 1. An optimal control sequence $u^{*}$ to (1) exists if and only if $\Omega\cap\mathcal{R}_{T}(x_{0})\neq\emptyset$.

Proof:

(Necessity) $\mathcal{R}_{T}(x_{0})$ includes all states reachable from $x_{0}$ at $t=T$ subject to the constraints in (1) (except $\Omega$), i.e., $x(T)\in\mathcal{R}_{T}(x_{0}),\forall u\in\Delta_{M}^{T}$. If $\Omega\cap\mathcal{R}_{T}(x_{0})=\emptyset$, no feasible solutions exist for Problem 1, i.e., Problem 1 is infeasible. Thus, the necessity of $\Omega\cap\mathcal{R}_{T}(x_{0})\neq\emptyset$ is obvious.

(Sufficiency) Note that the solution space $\mathbb{U}$ of Problem 1 is of finite size, which contains at most $M^{T}$ candidate solutions. Besides, $\Omega\cap\mathcal{R}_{T}(x_{0})\neq\emptyset$ implies that at least one feasible solution $u\in\mathbb{U}$ exists which can steer the BCN from $x_{0}$ to a terminal state $x_{f}\in\Omega$ at time $T$ under constraints. Moreover, Assumption 1 ensures that $J_{\textrm{T}}(u)$ is bounded from below, $\forall u\in\mathbb{U}$. A straightforward exhaustive enumeration of $\mathbb{U}$ can yield the optimal solution $u^{*}$ satisfying $J_{\textrm{T}}(u^{*})=\min_{u\in\mathbb{U}}J_{\textrm{T}}(u)$. ∎

Remark 7

If there are no constraints in Problem 1, i.e., $C_{x}\equiv\Delta_{N},C_{u}(\cdot)\equiv\Delta_{M},\Omega\equiv\Delta_{N}$, an optimal control sequence always exists, which is widely studied in existing work, e.g., [23, 17]. Note additionally that the optimal control sequence $u^{*}$ may not be unique.

V-B Time-Expanded Fixed-Time State Transition Graph

Fig. 1 shows that each edge in the STG corresponds to a state transition of the BCN, whose weight indicates the transition cost. This fact motivates us to connect Problem 1 to the shortest-path (SP) problem on the STG. However, in contrast to the standard SP problem in graph theory, Problem 1 poses three substantial challenges. First, the number of time steps is fixed to $T$, that is, we want only $T$-edge paths. Second, the stage cost function $g$ is time-dependent, indicating that the edge weights may vary with time. Finally, there is an additional terminal cost given by $h_{\textrm{T}}$. Consequently, the classic SP algorithms can no longer be applied. To overcome these obstacles, we get inspiration from the space-time network in dynamic transportation network studies [34] and propose a new graph called the Time-Expanded fixed-Time State Transition Graph (TET-STG), which attaches timestamp to state transitions by stretching the STG along the time dimension. Besides, we introduce a pseudo-state $\delta_{N}^{0}$ to handle terminal states and their costs. A formal definition is given below.

Definition 5

Consider Problem 1. The TET-STG $G_{\text{tet}}=(V,E,x_{0},T)$ is a weighted directed graph constructed by:

• •

$V=\cup_{t=0}^{T+1}V_{t}$, where $V_{t}=\mathcal{R}_{t}(x_{0}),\forall t\in[0,T-1]$, $V_{T}=\mathcal{R}_{T}(x_{0})\cap\Omega$, and $V_{T+1}=\{\delta_{N}^{0}\}$.

• •

$E=\cup_{t=0}^{T}E_{t}$*, where $E_{t}=\{(\delta_{N}^{i},\delta_{N}^{j})|\delta_{N}^{i}\in V_{t},\delta_{N}^{j}\in V_{t+1}\cap\mathcal{R}_{1}(\delta_{N}^{i})\},\forall t\in[0,T-1]$, and *

$E_{T}=\{(\delta_{N}^{i},\delta_{N}^{0})|\delta_{N}^{i}\in V_{T}\}$.

Note that a state may appear at multiple time points, but they are treated as distinctive vertices from a graph perspective. Denote the vertex representing the state $\delta_{N}^{i}$ at time $t$ by $\delta_{N,t}^{i}$. The weight of the edge $(\delta_{N,t}^{i},\delta_{N,t+1}^{j})\in E_{t}$ is

[TABLE]

where $U^{ij}$ is given in (9), and the unique vertex at time $T+1$ refers to the pseudo-state $\delta_{N,T+1}^{0}$. Denote the control that achieves the weight (cost) in (14) by $u_{t}^{ij}$:

[TABLE]

Remark 8

We note a slight abuse of notations in the above: $\mathcal{R}_{t}(x_{0})$ includes time information implicitly, and thus $V_{i}\cap V_{j}=\emptyset,\forall i\neq j$. Besides, recall Eq. (9): among the possibly nonunique control inputs that enable a state transition, we pick definitely the one of lowest cost in (14) for optimal control purpose. The role of the pseudo-state $\delta_{N,T+1}^{0}$ is to incorporate the terminal cost into the graph at the pseudo-time $T+1$.

Despite its seemingly complex definition, the TET-STG can be built handily by acquiring $V_{0},E_{0},V_{1},\cdots,E_{T-1},V_{T}$ successively similar to the BFS in Algorithm 1. In practical implementation, the one-step transition between states need to computed only once: supposing there is a transition $(\delta_{N}^{i},\delta_{N}^{j})$ in the STG, i.e., $\delta_{N}^{j}\in\mathcal{R}_{1}(\delta_{N}^{i})$, if we have a vertex $\delta_{N,t}^{i}\in V_{t}$, then there exists a succeeding vertex $\delta_{N,t}^{j}\in V_{t+1}$ and an edge $(\delta_{N,t}^{i},\delta_{N,t+1}^{j})\in E_{t}$ in the TET-STG.

Example 2

Consider the BCN $\Sigma_{1}$ (12) in Example 1 again. In addition to the constraints (13), Problem 1 is set up by $T=4,x_{0}=\{\delta_{8}^{1}\}$, and $\Omega=\{\delta_{8}^{2},\delta_{8}^{6}\}$. The costs are:

[TABLE]

where $Q=\textrm{diag}(2,3,1,0)$ and $R=\textrm{diag}(3,5,4,0,1,3,6,0)$.

Take the transition $(\delta_{8}^{3},\delta_{8}^{7})$ as an example. We have $U^{37}=\{\delta_{4}^{1},\delta_{4}^{3}\}$, which justifies the use of (15). The obtained TET-STG is shown in Fig. 2. Now it is clear that, though the stage cost function $g$ itself is time-dependent, the weight of each edge in the TET-STG becomes time-invariant after we expand the STG along the time axis. ∎

Recall Proposition 1. From the construction of the TET-STG, it is clear that an optimal control sequence exists if $V_{T}\neq\emptyset$, because each path from $\delta_{N,0}^{i_{0}}$ to $\delta_{N,T+1}^{0}$ yields a feasible solution to Problem 1, where $\delta_{N}^{i_{0}}=x_{0}$.

As aforementioned, we will transform Problem 1 into an SP problem in the TET-STG. First, note that each one-step transition in the TET-STG is already attained with a minimum cost thanks to (14) and (15). A direct consequence is the following lemma whose correctness is intuitive.

Lemma 3

Consider Problem 1 under Assumption 1. For any feasible control sequence $u=(\delta_{M}^{k_{0}},\delta_{M}^{k_{1}},\cdots,\delta_{M}^{k_{T-1}})$ that steers the BCN along a state trajectory $s=(\delta_{N}^{i_{0}},\delta_{N}^{i_{1}},\cdots,\delta_{N}^{i_{T}})$ with $x_{0}=\delta_{N}^{i_{0}}$ and $\delta_{N}^{i_{T}}\in\Omega$, there holds $J_{\textrm{T}}(u)\geq w(p)$, where $p$ is the corresponding path in the associated TET-STG,

[TABLE]

and $J_{\textrm{T}}(u)=w(p)$ holds if $\delta_{M}^{k_{t}}=u_{t}^{i_{t}i_{t+1}},\forall t\in[0,T-1]$.

Proof:

According to edge weights in (14) and the optimal one-step control in (15), the weight of the path $p$ in (16) is

[TABLE]

From the definition of $J_{\textrm{T}}$ in (1), we have

[TABLE]

We can see that the correctness of Lemma 3 is obvious. ∎

The following theorem establishes the connection between fixed-time optimal control and the SP in the TET-STG.

Theorem 1

Consider Problem 1 under Assumption 1 and suppose it is feasible. Given $x_{0}=\delta_{N}^{i_{0}}$, if an SP from $\delta_{N,0}^{i_{0}}$ to $\delta_{N,T+1}^{0}$ in the TET-STG $G_{\text{tet}}=(V,E,x_{0},T)$ is

[TABLE]

then the minimum value of the cost function (1) is $J_{\textrm{T}}^{*}=w(p^{*})$, and an optimal control sequence is

[TABLE]

Proof:

Suppose the solution space of Problem 1 is $\mathbb{U}\neq\emptyset$, and the set of paths from $\delta_{N,0}^{i_{0}}$ to $\delta_{N,T+1}^{0}$ in $G_{\text{tet}}$ is $\mathbb{P}$. There holds $\mathbb{P}\neq\emptyset$ because each state trajectory driven by $u\in\mathbb{U}$ corresponds to a path in $\mathbb{P}$. Moreover, we have $|\mathbb{P}|\leq|\mathbb{U}|$ since different control sequences may lead to identical state trajectories (see (9)). As we have shown in the proof of Proposition 2, $|\mathbb{U}|$ is finite, which implies $|\mathbb{P}|$ is also finite. Additionally, Assumption 1 and Eq. (14) ensures that all edge weights in $G_{\text{tet}}$ are bounded from below, which guarantees $w(p)$ is bounded from below for any $p\in\mathbb{P}$ because $p$ has exactly $T+1$ edges. Thus, there must exist a shortest path $p^{*}$.

Given any $u\in\mathbb{U}$, Lemma 3 tells that there exists $p\in\mathbb{P}$ such that $J_{\textrm{T}}(u)\geq w(p)\geq w(p^{*})$. Furthermore, for the control sequence $u^{*}\in\mathbb{U}$ in (18) that attains $p^{*}$, it leads to a state trajectory $s^{*}=(\delta_{N}^{i_{0}},\delta_{N}^{i_{1}^{*}},\cdots,\delta_{N}^{i_{T}^{*}})$ with $\delta_{N}^{i_{T}^{*}}\in\Omega$ according to the TET-STG in Definition 5. Thus, $u^{*}$ is a feasible solution to Problem 1. Furthermore, Lemma 3 states $J_{\textrm{T}}(u^{*})=w(p^{*})$. Hence, we have $J_{\textrm{T}}(u)\geq J_{\textrm{T}}(u^{*})$, $\forall u\in\mathbb{U}$. ∎

V-C Dynamic Programming (DP) in TET-STG

The problem following Theorem 1 immediately is how to locate an SP from $\delta_{N,0}^{i_{0}}$ to $\delta_{N,T+1}^{0}$ in the TET-STG efficiently (see (17)). Since the layered structure of the TET-STG doesn’t contain any cycles, the classic SP algorithms [33] such as the Floyd-Warshall algorithm and the Dijkstra’s algorithm can be applied directly. However, in view of the fixed horizon length in Problem 1, we propose a custom method based on dynamic programming (DP) to achieve better time efficiency, which can even beat the state-of-the-art Dijkstra’s algorithm (Remark 10).

The intuition behind our DP approach is that any sub-path of an SP is itself an SP as well [33]. Such optimal sub-structure is a strong indicator that DP based methods are applicable. The following theorem formalizes this idea.

Theorem 2

Consider Problem 1 under Assumption 1 and suppose it is feasible. In its associated TET-STG $G_{\text{tet}}=(V,E,x_{0},T)$ with $x_{0}=\delta_{N}^{i_{0}}$, let $F(\delta_{N,t}^{j})$ denote the weight of an SP from vertex $\delta_{N,0}^{i_{0}}$ to vertex $\delta_{N,t}^{j}$, $\forall t\in[0,T+1]$, and let $P(\delta_{N,t}^{j})$ denote the predecessors of $\delta_{N,t}^{j}$ in $G_{\text{tet}}$:

[TABLE]

Then $F(\delta_{N,t}^{j})$ can be obtained by the following recursion:

[TABLE]

for $t\in[1,T+1]$, and the base condition is $F(\delta_{N,0}^{i_{0}})=0$. If $p^{*}$ is an SP from $\delta_{N,0}^{i_{0}}$ to $\delta_{N,T+1}^{0}$, we have $w(p^{*})=F(\delta^{0}_{N,T+1})$.

Proof:

We first show the correctness of the recursion (20) by induction. Due to the layered structure of $G_{\text{tet}}$, any vertex $\delta_{N,t}^{j}\in V$ can only be reached from a certain vertex $\delta_{N,t-1}^{i}\in V$ in one step, $t\geq 1$. Assume $F(\delta_{N,t-1}^{i}),t\geq 1,$ represents the minimum weight from vertex $\delta_{N,0}^{i_{0}}$ to vertex $\delta_{N,t-1}^{i}$. For any path from $\delta_{N,0}^{i_{0}}$ to $\delta_{N,t}^{j}$ that passes $\delta_{N,t-1}^{i}$, its minimum weight is obviously $F(\delta_{N,t-1}^{i})$ + $w(\delta_{N,t-1}^{i},\delta_{N,t}^{j})$. Eq. (20) examines all such predecessors of $\delta_{N,t}^{j}$, and the minimum value $F(\delta_{N,t}^{j})$ is clearly the minimum weight of any path from $\delta_{N,0}^{i_{0}}$ to $\delta_{N,t}^{j}$. Besides, for $t=0$, the base case $F(\delta_{N,0}^{i_{0}})=0$ is clearly true. Thus, we have verified that $F(\delta^{j}_{N,t})$ is the minimum weight of any path from $\delta_{N,0}^{i_{0}}$ to $\delta^{j}_{N,t}$. It implies directly that $F(\delta^{0}_{N,T+1})$ is the weight of the SP $p*$ from $\delta_{N,0}^{i_{0}}$ to $\delta_{N,T+1}^{0}$. ∎

Remark 9

Eq. (20) is essentially a form of the Bellman optimality equation, a widely applied tool in optimal control [3], though we state it from a graph-theoretical perspective.

Combing Theorem 1 and 2, we get the minimum cost $J_{\textrm{T}}^{*}$ for Problem 1 by $J_{\textrm{T}}^{*}=w(p^{*})=F(\delta^{0}_{N,T+1})$. Nonetheless, we are more interested in the optimal control sequence $u^{*}$ that attains $J_{\textrm{T}}^{*}$. The key is to record the minimizer to (20): if $\delta_{N,t}^{j}$ is in an SP, then the minimizer $\delta_{N,t-1}^{i^{*}}$ must be in the SP as well. At last, $p^{*}$ can be reconstructed accordingly, followed by $u^{*}$ acquired with (18). Algorithm 2 implements the DP method by Theorem 2 to solve $F(\delta_{N,t}^{j})$ and to reconstruct the optimal control sequence. Note that, for maximal time efficiency, the predecessors (19) of each vertex are stored while building $G_{\text{tet}}$.

Time complexity analysis of Algorithm 2. $G_{\text{tet}}$ has $T+2$ layers with at most $Z\vcentcolon=|\mathcal{R}(x_{0})|$ vertices per layer and at most $MZ$ edges between any two successive layers. Like Algorithm 1, $G_{\text{tet}}$ is constructed via BFS in linear time $O(TMZ)$. The core function of Algorithm 2, ShortestPath, implements top-down DP via the memoization technique [33], i.e., storing results into $F$ and retrieving the cached results if same inputs recur. Memoization can avoid repetitive computation [33], and accordingly each edge of $G_{\text{tet}}$ is processed only once. ShortestPath in Line 19 thus takes time $O(TMZ)$. The remaining construction of $u^{*}$ runs in $O(T)$. The overall worst-case time complexity of Algorithm 2 is hence $O(TMZ)$, or equivalently, $O(TMN)$, because we always have $Z\leq N$.

Remark 10

The recent work [23] establishes a weighted directed graph to formulate the $k$-edge shortest path problem as well. However, neither the time-dependent stage cost function nor the terminal constraint set is considered in [23]. Moreover, the work [23] directly applies the standard Dijkstra’s algorithm, which is designed for general SP problems and less efficient than our Algorithm 2 in solving Problem 1, whose running time is $O(TMN+TN\log(TN))$ instead.

Example 3

Recall Example 2. Running Algorithm 2, we get the minimum value of the cost function $J_{\textrm{T}}^{*}=11$, and the optimal control sequence $u^{*}=(\delta_{4}^{4},\delta_{4}^{3},\delta_{4}^{4},\delta_{4}^{3})$. The correctness of this result can be easily verified by inspecting Fig. 2 and enumerating all paths from $\delta_{4,0}^{1}$ to $\delta_{4,5}^{0}$. ∎

VI Solve Problem 2 via Dijkstra’s Algorithm

In this section, the existence of optimal solutions to Problem 2 is first examined. Then, we divide Problem 2 into two cases depending on whether cost functions are time-dependent. Both cases will be conquered by Dijkstra’s algorithm, but different graph structures are constructed to maximize efficiency.

VI-A Existence of Optimal Solutions

Like Problem 1, we use the terminal costs and the stage costs of Problem 2 as the edge weights of specific graphs. However, the major difference is that the number of state transitions in Problem 2 is not fixed. Consequently, the condition that no negative-weight cycles exist in any state trajectory from $x_{0}$ to $x_{d}\in\Omega$ is mandatory [33]; otherwise, the cost $J$ can always be reduced by repeating a negative-weight cycle, and no SP exists. We thus require the following conditions to guarantee existence of optimal solutions to Problem 2.

Assumption 2

The cost functions in Problem 2 satisfy three conditions: (i) $h$ is bounded from below; (ii) $g$ is nonnegative; (iii) $h$ and $g$ are both nondecreasing with respect to time $t$, i.e., $h(\delta_{N}^{i},t_{2})\geq h(\delta_{N}^{i},t_{1}),\forall t_{2}>t_{1},\forall\delta_{N}^{i}\in\Omega$, and $g(\delta_{N}^{i},\delta_{M}^{k},t_{2})\geq g(\delta_{N}^{i},\delta_{M}^{k},t_{1}),\forall t_{2}>t_{1},\forall(\delta_{N}^{i},\delta_{M}^{k})\in\Delta_{N}\times\Delta_{M}$.

The rationality of condition (i) is obvious, just like Problem 1, to ensure a finite optimal value. Condition (ii) can be technically relaxed to the nonexistence of negative-weight cycles, though it would be quite difficult to verify such a condition in practice because edge weights (i.e., stage costs) can vary with time. We can justify condition (iii) intuitively by imaging a special scenario. Suppose a state trajectory from $x_{0}$ to $x_{d}\in\Omega$ contains a cycle of zero weight. Then a possible result is that the more cycling the BCN does along this cycle, the more the cost criterion $J$ can be reduced, once $g$ or $h$ can decrease as time passes. Consequently, the optimal control sequence does not have a finite length. Note that if $g$ and $h$ do not depend on $t$, which is the most common case in the literature, condition (iii) is satisfied naturally.

Remark 11

Following [17, 16], we can always assume $h$ is nonnegative without affecting the optimal solution. Supposing $h$ is bounded from below by a constant $B_{h}$, the optimal control sequence to (2) is the same one that minimizes $J^{\prime}(u)=J(u)-B_{h}$, and we have a nonnegative terminal cost function now: $h^{\prime}(x(K),K)=h(x(K),K)-B_{h}\geq 0$.

The following proposition confirms the existence of an optimal solution to Problem 2 under the above conditions.

Proposition 2

Consider Problem 2 under Assumption 2. There exists an optimal control sequence $u^{*}$ satisfying $|u^{*}|<|\mathcal{R}(x_{0})|$ that minimizes (2) if and only if $\Omega\cap\mathcal{R}(x_{0})\neq\emptyset$.

Proof:

(Necessity) Since $\mathcal{R}(x_{0})$ denotes all states reachable from $x_{0}$, no state $x_{d}\in\Omega$ can be reached if $\Omega\cap\mathcal{R}(x_{0})=\emptyset$.

(Sufficiency) If $\Omega\cap\mathcal{R}(x_{0})\neq\emptyset$, there must exist control sequences that steer the BCN from $x_{0}$ to a state $x_{d}\in\Omega$. Suppose such a control sequence is u=\big{(}u(t)\big{)}_{t=0}^{k-1}, and the resultant trajectory is s=\big{(}x(t)\big{)}_{t=0}^{k}, where $x(0)=x_{0}$ and $x(k)=x_{d}$. Furthermore, we claim that if $|u|=k\geq|\mathcal{R}(x_{0})|$, there must exist a shorter control sequence $\bar{u}$ such that $|\bar{u}|<|\mathcal{R}(x_{0})|$ and $J(\bar{u})\leq J(u)$. This claim is proved below.

$|u|\geq|\mathcal{R}(x_{0})|$ implies that $|s|>|\mathcal{R}(x_{0})|$, which means that $s$ contains repetitive states because $x\in\mathcal{R}(x_{0})$ for any state $x$ in $s$. Assume one such pair of repetitive states is $x(i)=x(j)$, $0\leq i<j\leq k$, and thus $c=(x(i),x(i+1),\cdots,x(j))$ is a circular sub-trajectory. We can remove this cycle (except $x(j)$) from $s$, and obviously the remaining states $s^{\prime}=(x(0),\cdots,x(i-1),x(j),x(j+1),\cdots,x(k-1),x(k))$ still constitute a trajectory from $x_{0}$ to $x_{d}$ driven by a shortened control sequence $u^{\prime}=(u(0),\cdots,u(i-1),u(j),u(j+1),\cdots,u(k-1))$. Let $r\vcentcolon=j-i>0$, and it follows that

[TABLE]

Condition (ii) and (iii) in Assumption 2 guarantee (VI-A) is nonnegative, i.e., $J(u^{\prime})\leq J(u)$. The above cycle elimination procedure can be repeated until a control sequence $\bar{u}$ satisfying $|\bar{u}|<|\mathcal{R}(x_{0})|$ is obtained. We have certainly $J(\bar{u})\leq J(u)$.

The above claim implies that it is enough to consider the candidate set $\bar{\mathbb{U}}=\{u||u|<|\mathcal{R}(x_{0})|,u\in\mathbb{U}\}$ for an optimal solution, where $\mathbb{U}$ is the feasible set of Problem 2. Obviously, the set $\bar{\mathbb{U}}$ is finite, and Assumption 2 guarantees $J(u)$ is bounded from below, $\forall u\in\bar{\mathbb{U}}$. Therefore, an optimal solution $u^{*}\in\bar{\mathbb{U}}$ to Problem 2 must exist such that $|u^{*}|<|\mathcal{R}(x_{0})|$. ∎

VI-B Case 1: Time-Invariant Stage Cost and Terminal Cost

In this case, neither $g$ nor $h$ of (2) depends on time $t$: the STG becomes a static graph, whose edge weights are permanently fixed. Furthermore, if there is only one destination state with zero terminal cost, Problem 2 degrades to a standard SP problem in the STG. This simplest case has been solved in [14, 15, 23]. We address the more general problems here, where multiple destination states with non-zero terminal costs are allowed. Following the same idea in solving Problem 1, we introduce an extra pseudo-state $\delta_{N}^{0}$ as well as the terminal set into the STG and term the new graph STG+.

Definition 6

Consider Case 1 of Problem 2. The extended state transition graph STG+, denoted by $G^{+}=(V,E,x_{0})$, is an extension of the STG constructed as follows:

Build the STG $G=(V,E,x_{0})$ by Definition 4, and the weights are assigned like (14), though time-independent:

[TABLE] 2. 2.

Add into $G$ the pseudo-state $\delta_{N}^{0}$ and its incoming edges,

[TABLE]

where $E^{0}=\{(\delta_{N}^{i},\delta_{N}^{0})|\delta_{N}^{i}\in V\cap\Omega\}$ with weights:

[TABLE]

Here $\alpha\in\mathbb{N}$ is an arbitrary integer used as a placeholder.

Since all edge weights of $G^{+}$ are time-invariant, we denote the weight and the associated control input of each edge $(\delta_{N}^{i},\delta_{N}^{j})\in E$ by $w^{ij}$ and $u^{ij}$ respectively. Of course, the incoming edges of $\delta_{N}^{0}$ needs on control input.

Example 4

Consider Example 1 again: the BCN (12) is subject to the constraints (13). Suppose the desired terminal states are $\Omega=\{\delta_{8}^{3},\delta_{8}^{4}\}$, and the initial state is $x_{0}=\delta_{8}^{7}$. We assign the following time-invariant cost:

[TABLE]

where $Q=\textrm{diag}(2,5,1,4,1,3,6,0)$ and $R=\textrm{diag}(0,3,1,4)$. The STG+of this case can be easily built on basis of the STG in Fig. 1 following Definition 6, which is shown in Fig. 3. As we see, the STG+is akin to the STG but with an additional pseudo-state and edge weights assigned. ∎

To transform Problem 2 in this case into an SP problem in the STG+, we first give the following lemma.

Lemma 4

Consider Case 1 of Problem 2 under Assumption 2 with time-independent $g$ and $h$. Given any control sequence $u=(\delta_{M}^{k_{0}},\delta_{M}^{k_{1}},\cdots,\delta_{M}^{k_{K-1}})$ of an unknown length $K$ that steers the BCN along a state trajectory $s=(\delta_{N}^{i_{0}},\delta_{N}^{i_{1}},\cdots,\delta_{N}^{i_{K}})$ with $x_{0}=\delta_{N}^{i_{0}}$ and $\delta_{N}^{i_{K}}\in\Omega$, there holds $J(u)\geq w(p)$, where $p$ is the corresponding path in the associated STG+,

[TABLE]

and $J(u)=w(p)$ holds if $\delta_{M}^{k_{t}}=u^{i_{t}i_{t+1}},\forall t\in[0,K-1]$.

Proof:

With the edge weights given by Definition 6, the weight of the path $p$ is

[TABLE]

Note that $u^{i_{t}i_{t+1}}$ is the optimal control input to transit the BCN from $\delta_{N}^{i_{t}}$ to $\delta_{N}^{i_{t+1}}$ in one step, i.e., $g(\delta_{N}^{i_{t}},u^{i_{t}i_{t+1}},t)\leq g(\delta_{N}^{i_{t}},\delta_{M}^{k_{t}},t),\forall t\in[0,K-1]$. Comparing $w(p)$ with the cost $J(u)$ in (2), Lemma 4 is obviously true. ∎

The following theorem relates Case 1 of Problem 2 to the SP problem in an STG+.

Theorem 3

Consider Case 1 of Problem 2 under Assumption 2 with time-independent $g$ and $h$, and suppose it is feasible. In the STG+ $G^{+}=(V,E,x_{0}),x_{0}=\delta_{N}^{i_{0}}$, there exists an SP p^{*}=\big{<}\delta_{N}^{i_{0}},\delta_{N}^{i_{1}^{*}},\delta_{N}^{i_{2}^{*}},\cdots,\delta_{N}^{i_{K-1}^{*}},\delta_{N}^{i_{K}^{*}},\delta_{N}^{0}\big{>},K\leq|V|-2, such that the minimum value of the cost function (2) is $J^{*}=w(p^{*})$. The corresponding optimal control sequence is $u^{*}=\{u^{i_{0}i_{1}^{*}},u^{i_{1}^{*}i_{2}^{*}},\cdots,u^{i_{K-1}^{*}i_{K}^{*}}\}$.

Proof:

Since the problem is feasible, there exists control sequences that steers the BCN from the initial state $x_{0}$ to a destination state $x_{d}\in\Omega$, implying that there exist paths from $\delta_{N}^{i_{0}}$ to $\delta_{N}^{0}$ in $G^{+}$. Note that $G^{+}$ has no negative cycles because (i) Assumption 2 states $g$ is nonnegative and (ii) the pseudo-state $\delta_{N}^{0}$ has only incoming edges, whose weights assigned by $h$ may be negative though. With this fact, we can easily show that an SP $p^{*}$ from $\delta_{N}^{i_{0}}$ to $\delta_{N}^{0}$ with at most $|V|$ vertices exists following the cycle elimination procedure in Proposition 2. In fact, this is a fundamental theorem in graph theory [33].

Recall Lemma 4: given any feasible control sequence $u$, there exists a path $p$ from $\delta_{N}^{i_{0}}$ to $\delta_{N}^{0}$ in $G^{+}$ such that $J(u)\geq w(p)\geq w(p^{*})=J(u^{*})$. The above theorem holds clearly. ∎

Theorem 3 has reduced Problem 2 with time-invariant costs to a regular single-pair SP problem [33] from $x_{0}$ to $\delta_{N}^{0}$ in the STG+. Since $g$ is nonnegative and the pseudo-state $\delta_{N}^{0}$ has only incoming edges, We have discussed in Remark 11 that we can always assume $h$ is nonnegative without loss of generality, though it is not mandatory. Since $g$ is also nonnegative, the STG+has only nonnegative edge weights. The fastest known SP algorithm for such graphs is Dijkstra’s algorithm [33].

We make two modifications to the normal implementation of Dijkstra’s algorithm for this optimal control problem. First, like Algorithm 2, the vertices that compose the SP are recorded to reconstruct the optimal control sequence later. Second, we terminate the search process once vertex $\delta_{N}^{0}$ is reached because we are only interested in the SP from $x_{0}$ to $\delta_{N}^{0}$. Algorithm 3 presents the modified Dijkstra’s algorithm to solve Case 1 of Problem 2. Since Dijkstra’s algorithm is a well-established SP algorithm in graph theory (see [33, Chapter 24] for details), the proof of its correctness is omitted here.

Time complexity analysis of Algorithm 3. A key data structure in Dijkstra’s algorithm is the priority queue [33], in which each item has a priority and the one with highest (or lowest) priority is first served. If the priority queue is implemented with a Fibonacci heap, Dijkstra’s algorithm has a running time of $O(|E|+|V|\log|V|)$ [33]. In the beginning, just like the STG, the construction of the STG+takes time $O(M|V|)$. The more complicated Dijkstra’s part (Line 2 to 16) runs in $O(M|V|+|V|\log|V|)$ accordingly. Finally, since the SP found by Dijkstra’s algorithm contains at most $|V|$ vertices, the construction of $u^{*}$ (Line 17 to 22) runs in $O(|V|)$. Overall, the worst-case time complexity of Algorithm 3 is $O(M|V|+|V|\log|V|)$, or equivalently, $O(MN+N\log N)$.

Example 5

We test Algorithm 3 with Example 4. The optimal value is $J^{*}=13$ and the optimal control sequence is $u^{*}=(\delta_{4}^{1},\delta_{4}^{3},\delta_{4}^{1})$. The state trajectory of the BCN is thus $s^{*}=(\delta_{8}^{7},\delta_{8}^{5},\delta_{8}^{2},\delta_{8}^{4})$. It is easy to verify in Fig. 3 that an SP from $\delta_{8}^{7}$ to $\delta_{8}^{0}$ is p^{*}=\big{<}\delta_{8}^{7},\delta_{8}^{5},\delta_{8}^{2},\delta_{8}^{4},\delta_{8}^{0}\big{>} and $w(p^{*})=13$. ∎

VI-C Case 2: Time-Variant Stage Cost and Terminal Cost

As aforementioned, the classic SP algorithms will not work once the edge weights may vary with time. To the best of knowledge, there are still no published studies on fixed-destination optimal control of BCNs with time-varying costs. Recall the TET-STG proposed in Section V, and we naturally attempt to handle this time-variant case for Problem 2 in a similar way. However, one immediate difficulty is that, unlike Problem 1, the horizon length is not known beforehand, which prevents the reuse of the DP-based Algorithm 2 directly.

Hopefully, we may resort to Proposition 2 to overcome this obstacle: it is sufficient to consider only control sequences of size less than $|\mathcal{R}(x_{0})|$ to find the optimal one, though the exact length remains unknown. On the other hand, Algorithm 2 works once a horizon length $T$ is given. A straightforward workaround that reuses Algorithm 2 to solve this case thus comes to our mind as follows.

Procedure 1

Solve Case 2 by Reusing Algorithm 2.

• •

Step 1. Given $T\in[1,|\mathcal{R}(x_{0})|-1]$, transform this case into Problem 1 by setting $h_{\textrm{T}}(x(T))=h(x(T),T)$.

• •

Step 2. Solve the problem obtained above with Algorithm 2 to get an optimal control sequence $u_{T}^{*},|u_{T}^{*}|=T$.

• •

Step 3. Repeat Step 1 and Step 2 for all possible $T$’s, and finally yield $u^{*}=\operatorname*{arg\,min}_{u_{T}^{*}}(J(u_{T}^{*}))$.

The enumeration of all possible horizon lengths like Procedure 1 is essentially the idea underlying the algebraic approach [14] for minimum-energy control towards a given target state, though it only considers time-invariant stage cost. A similar idea is also adopted in [20] to detect the minimum average-weight cycle in the input-state space. However, such a somewhat brute-force method is still inevitably computationally expensive, e.g., the running time is $O(N^{4})$ in [14, 20]. In our case, even though each subproblem for a specific $T$ can be solved by the more efficient Algorithm 2, the overall time complexity of Procedure 1 is still as high as $O(MN^{3})$.

Since we only care about control sequences that has a size less than $|\mathcal{R}(x_{0})|$, to further reduce the computational burden, we devise another approach by adapting the TET-STG to an unknown but limited horizon length. More interestingly, Algorithm 3 initially developed for Case 1 can be reused on the resultant graph. We call this new graph a Time-Expanded fixed-Destination State Transition Graph (TED-STG). It construction is similar to the TET-STG detailed as follows.

Definition 7

Consider Case 2 of Problem 2. The TED-STG $G_{\text{ted}}=(V,E,x_{0})$ is a weighted directed graph constructed by:

• •

$V=\cup_{t=0}^{Z}V_{t}$, where $V_{t}=\mathcal{R}_{t}(x_{0}),\forall t\in[0,Z-1]$, $V_{Z}=\{\delta_{N}^{0}\}$, and $Z\vcentcolon=|\mathcal{R}(x_{0})|$ is the reachable set size.

• •

$E=\cup_{t=0}^{Z-2}E_{t}\cup E^{0}$, where $E_{t}=\{(\delta_{N}^{i},\delta_{N}^{j})|\delta_{N}^{i}\in V_{t},\delta_{N}^{j}\in V_{t+1}\cap\mathcal{R}_{1}(\delta_{N}^{i})\},\forall t\in[0,Z-2]$, and connect the terminal states in each layer to the pseudo-state $\delta_{N}^{0}$ by $E^{0}=\{(\delta_{N,t}^{i},\delta_{N}^{0})|\delta_{N}^{i}\in\Omega\cap V_{t},t\in[0,Z-1]\}$.

A vertex $\delta_{N,t}^{i}$ above refers to the state $\delta_{N}^{i}$ at time $t$. The weight of the edge $(\delta_{N,t}^{i},\delta_{N,t+1}^{j})\in E_{t},\forall t\in[0,Z-2],$ is

[TABLE]

and the control enabling the transition from $\delta_{N}^{i}$ to $\delta_{N}^{j}$ at time $t$ is also $u_{t}^{ij}$ in (15). The weight of edges in $E^{0}$ is

[TABLE]

Remark 12

Since the destination (terminal) states may be reached at any time in an optimal trajectory, we package all such possibilities into $E^{0}$. This greatly simplifies the problem: we only need to find an optimal path from $\delta_{N,0}^{i_{0}}$ to $\delta_{N}^{0}$. Accordingly, the time to reach the pseudo-state $\delta_{N}^{0}$ is not known in advance, and that is why it has no time subscript.

Example 6

We use the BCN (12) in Example 1 to illustrate the TED-STG. Problem 2 is set up by $x_{0}=\{\delta_{8}^{1}\}$ and $\Omega=\{\delta_{8}^{6}\}$. The time-variant stage cost and terminal cost are:

[TABLE]

where $Q=\textrm{diag}(2,3,1,5+t)$ and $R=\textrm{diag}(3,2t,4,0,1,5+t,6,0)$. The TED-STG of this task is shown in Fig. 4. Note that the size of $\mathcal{R}(x_{0})$ is $Z=7$ here. We consider only one destination state to facilitate illustration. ∎

Just like the TET-STG, the edge weights in a TED-STG do not change with time, though the cost functions $g$ and $h$ are themselves time-dependent. The cost of a state trajectory is related with the weight of a path in the TED-STG as follows.

Lemma 5

Consider Case 2 of Problem 2 under Assumption 2. Given any control sequence $u=(\delta_{M}^{k_{0}},\delta_{M}^{k_{1}},\cdots,\delta_{M}^{k_{K-1}})$ of length $K<|\mathcal{R}(x_{0})|$ that steers the BCN along a state trajectory $s=(\delta_{N}^{i_{0}},\delta_{N}^{i_{1}},\cdots,\delta_{N}^{i_{K}})$ with $x_{0}=\delta_{N}^{i_{0}}$ and $\delta_{N}^{i_{K}}\in\Omega.$ There holds $J(u)\geq w(p)$, where $p$ is the corresponding path from $\delta_{N,0}^{i_{0}}$ to $\delta_{N}^{0}$ in the associated TED-STG,

[TABLE]

and $J(u)=w(p)$ holds if $\delta_{M}^{k_{t}}=u_{t}^{i_{t}i_{t+1}},\forall t\in[0,K-1]$.

Lemma 5 can be proved easily in almost the same way as Lemma 3, which is omitted here. Now we are ready to solve Case 2 of Problem 2 by converting it to an SP problem in the TED-STG through the following theorem.

Theorem 4

Consider Case 2 of Problem 2 under Assumption 2 with time-dependent $g$ and $h$, and suppose it is feasible. Given the initial state $x_{0}=\delta_{N}^{i_{0}}$, there exists an SP p^{*}=\big{<}\delta_{N,0}^{i_{0}},\delta_{N,1}^{i_{1}^{*}},\cdots,\delta_{N,K-1}^{i_{K-1}^{*}},\delta_{N,K}^{i_{K}^{*}},\delta_{N}^{0}\big{>},K<|\mathcal{R}(x_{0})|, in the TED-STG $G_{\text{ted}}=(V,E,x_{0})$ such that the minimum value of the cost function (2) is $J^{*}=w(p^{*})$. The corresponding optimal control sequence is $u^{*}=\{u_{t}^{i_{0}i_{1}^{*}},u_{t}^{i_{1}^{*}i_{2}^{*}},\cdots,u_{t}^{i_{K-1}^{*}i_{K}^{*}}\}$.

Proof:

Since the problem is feasible, Proposition 2 implies there exists an optimal control sequence shorter than $|\mathcal{R}(x_{0})|$. Thus, we can search $\mathbb{U}^{\prime}=\{u||u|<|\mathcal{R}(x_{0})|,u\in\mathbb{U}\}$ for an optimal one, where $\mathbb{U}$ is the feasible set of Problem 2. Now recall Lemma 5: given any $u\in\mathbb{U}^{\prime}$, there exists a path $p$ from $\delta_{N,0}^{i_{0}}$ to $\delta_{N}^{0}$ in the TED-STG such that $J(u)\geq w(p)\geq w(p^{*})$. Besides, with the weight given by Definition 7, it is obvious that $w(p^{*})=J(u^{*})$. Thus, the proof is complete. ∎

Theorem 4 transforms Problem 2 with time-dependent cost functions into a standard single-pair (i.e., $\delta_{N,0}^{i_{0}}$ to $\delta_{N}^{0}$) SP problem in the TED-STG. As discussed in Remark 11, we assume that $h\geq 0$ without loss of generality. Since all edges of this graph have nonnegative weights, we can apply Dijkstra’s algorithm again to identify the SP, which is the same as Algorithm 3 except that the TED-STG is used instead of the STG+. We detail this algorithm in the online supplementary material on ArXiv111Refer to https://arxiv.org/abs/1908.02019 for the supplementary material. (Algorithm LABEL:S-alg:_dijkstra_case2) to conserve space here.

Time complexity analysis of Algorithm LABEL:S-alg:_dijkstra_case2. Definition 7 tells that there are $Z\vcentcolon=|\mathcal{R}(x_{0})|$ layers in the TED-STG. Each layer has at most $Z$ vertices, and at most $MZ$ edges exist between two consecutive layers. Besides, there are typically only few destination states, i.e., the number of incoming edges of the pseudo-state $\delta_{N}^{0}$ is $O(Z)$. We thus have $|V|=O(Z^{2})$ and $|E|=O(MZ^{2})$. Like a TET-STG, the TED-STG can be built quickly in linear time, i.e., $O(|V|+|E|)=O(MZ^{2})$. Thus, the time complexity of Algorithm LABEL:S-alg:_dijkstra_case2 is dominated by the Dijkstra’s SP part. To conclude, Algorithm LABEL:S-alg:_dijkstra_case2 runs in time $O(MZ^{2}+Z^{2}\log Z^{2})=O(Z^{2}(M+2\log Z))$, which is much faster than the naive Procedure 1. Since we have $Z\leq N$, the time complexity is equivalent to $O(N^{2}(M+2\log N))$.

Example 7

Recall Example 6 and its TED-STG in Fig. 4. Algorithm LABEL:S-alg:_dijkstra_case2 yields the following results: $J^{*}=14$ and $u^{*}=(\delta_{4}^{3},\delta_{4}^{3},\delta_{4}^{3},\delta_{4}^{1})$. The corresponding SP in the TED-STG is $p^{*}=(\delta_{8,0}^{1},\delta_{8,1}^{4},\delta_{8,2}^{7},\delta_{8,3}^{5},\delta_{8,4}^{6},\delta_{8}^{0})$. ∎

VII Comparison with Related Work

As we have reviewed in Section I, unlike our algorithms which target the most general problems, most existing methods are developed for certain special cases of Problem 1 or Problem 2. We therefore categorize various optimal control tasks according to their characteristics to facilitate comparison. Their time complexity is summarized in Table I, where, as always, $N=2^{n}$ and $M=2^{m}$ for a $n$-state, $m$-input BCN, and $T$ denotes the fixed horizon length in Problem 1.

To better interpret Table I, note that we can always assume $M\leq N$ because a state can transit to at most $N$ succeeding states regardless of the number of control inputs. In fact, we usually have $m<n$ and thus $M\ll N$ in practice especially for large networks. For example, it can be enough to steer the whole network by controlling only a fraction of the nodes [30, 35]. In short, Table I shows that our graph-theoretical approach can accomplish higher time efficiency than most existing approaches, and only methods in [16] and [23] share the same time complexity as ours for particular problems. Notably, if there are time-variant costs in Problem 2, only Algorithm LABEL:S-alg:_dijkstra_case2 can handle it to the best of our knowledge. In summary, though we target FHOC problems ambitiously in their most general form, the computational efficiency of our approach is still superior to that of most existing work.

Note that the time complexity listed in Table I refers to the worst-case one, which indicates the longest running time of an algorithm given any possible input. By convention, the worst-case running time is used to measure the efficiency of algorithms [33, 17, 22]. A noteworthy point is that all the algebraic approaches in Table I, i.e., all existent work except [23], have their average-case time complexity equal to the worst-case one, because they essentially operate on matrices of identical sizes irrespective of the initial state $x_{0}$ and the size of its reachable set $\mathcal{R}(x_{0})$. By contrast, as shown in time complexity analysis of our algorithms, the actual size of the graph depends on $\mathcal{R}(x_{0})$, which is typically a small subset of the state space, while the algebraic methods always consider all the $N$ states. Additionally, in constraint handling, our approach excludes the undesirable states and transitions completely from the graph, but most algebraic approaches simply assign them infinitely large cost values and still involve them in subsequent operations. Consequently, our graph-theoretical approach attains potentially lower average-case time complexity than algebraic methods in practice, like [16], even though they share the same worst-case complexity.

VIII A Biological Benchmark Example

In this section, we focus on a larger BCN, i.e., the Ara operon gene regulatory network in E. coli that is responsible for sugar metabolism regulation [36, 22]. The main purpose is to compare the computational efficiency of various approaches. The Boolean model of this network is given in Table II, where the target nodes indicate the 9 state variables, and the 4 control inputs are $\{A_{e},A_{em},A_{ra_{-}},G_{e}\}$. Interested readers may consult [36] for the biological meaning of these variables. Using the STP, we can get the ASSR of this network with a network transition matrix $L\in\mathcal{L}_{N\times MN},N=512,M=16$.

To facilitate comparison with existing studies, we consider two tasks widely studied in the literature, i.e., the minimum-energy control and the minimum-time control. In both tasks, no constraints are enforced because only few existing methods consider state or input constraints. Suppose the initial state is $x_{0}=\delta_{512}^{9}$ and the desired state is $x_{d}=\delta_{512}^{410}$ for both tasks. All algorithms were implemented with Python 3.7. We did experiments on a desktop PC equipped with a 3.4 GHz Core i7-3770 CPU, 16 GB RAM, and Windows 10.

VIII-A Task 1: Minimum-Energy Control

In this task, we aim to transfer the BCN from the initial state $x_{0}$ to the desired state $x_{d}$ at a prespecified time point with least energy consumption [14]. This task is presented as an instance of Problem 1 with $T=10$ and $\Omega=\{x_{d}\}$. Since most methods can only deal with time-invariant costs, we use a stage cost function like that in [22] to evaluate the energy consumed by each state transition, and set zero terminal costs:

[TABLE]

where $X(t)=[x_{1}(t),x_{2}(t),\cdots,x_{9}(t)]^{\top}$ and $U(t)=[u_{1}(t),u_{2}(t),u_{3}(t),u_{4}(t)]^{\top}$. The two weight vectors are $\mathcal{A}=[0,16,40,44,28,28,28,48,44]$, and $\mathcal{B}=[0,48,28,48]$.

Running Algorithm 2, we get the following results.

• •

The reachable set of $x_{0}$, i.e., $\mathcal{R}(x_{0})$, has only 108 states, though the complete state state has 512 states in total.

• •

The minimum value of $J_{\textrm{T}}(u)$ in (1) is $J_{\textrm{T}}^{*}=1108$.

• •

The optimal control sequence is $u^{*}=\delta_{16}(16,16,16,16,16,16,8,5,6,14)$, and the resultant state trajectory is $s^{*}=\delta_{512}(9,457,463,480,480,480,480,352,312,288,410).$

We illustrate $\mathcal{R}(x_{0})$ and highlight the above optimal state trajectory $s^{*}$ in Fig. 5. We tested the other methods, and they all obtained the same optimal value $J_{T}^{*}$. The running time of different methods is listed in Table III. As we have expected, the algorithm in [16] takes more time than ours in practice, though they have identical worst-case time complexity. The main reason is that the former always evaluates the whole state space, while our algorithm only focuses on the reachable states $\mathcal{R}(x_{0})$. Another method [23] shares the same view as ours, but it depends on the slower Dijkstra’s algorithm rather than the more efficient DP-based Algorithm 2. Overall, our approach attains the shortest running time to complete Task 1.

VIII-B Task 2: Minimum-Time Control

This task is also referred to as time-optimal control, that is, to steer the BCN from $x_{0}$ to a desired state $x_{d}$ as fast as possible [15]. It is easy to specialize Problem 2 for this task: just let $g(\cdot,\cdot,\cdot)\equiv 1$, i.e., unit time for each state transition, $h(\cdot,\cdot)\equiv 0$, and $\Omega=\{x_{d}\}$. In this setting, the optimal value $J^{*}$ of (2) is obviously the minimum time taken from $x_{0}$ to $x_{d}$. Since $g$ and $h$ are both time-independent, we can tackle this task with Algorithm 3. The results are given below.

• •

The optimal value is $J^{*}=3$, i.e., the BCN is transferred from $x_{0}$ to $x_{d}$ in at least 3 steps.

• •

The optimal control sequence is $u^{*}=\delta_{16}(1,2,14)$, which leads to the state trajectory $s^{*}=\delta_{512}(9,41,15,410)$.

Like Task 1, we sketch the above minimum-time state trajectory in Fig. 5. The other methods yielded the same minimum time, and their running time is listed in Table III. In this task, the method in [23] is essentially identical to ours, both depending on Dijkstra’s algorithm to find an SP, and they are the fastest ones, taking far less time than the other two. Note that we have optimized it by building the STG efficiently with Algorithm 1 instead of the expensive algebraic method originally used in [23]. The algorithm in [14] is extremely slow here mainly because it examines exhaustively all SP’s of length ranging from 1 to $N$ to find the shortest one.

Remark 13

The running time of some algebraic approaches, like [17] and [15], might be further reduced using advanced numerical routines, since they depend heavily on matrix operations, and the involved matrices are typically sparse. Nevertheless, the results in Table III still demonstrate the supreme efficiency of our approach with an advantage of several orders of magnitude. Besides, the method in [23] is closest to ours, but it can only solve a small subset of problems investigated in this study. The Python implementation of the proposed approach and existing algorithms is available at GitHub https://github.com/ShuhuaGao/FHOC.

IX Conclusion

This paper focused on FHOC of BCNs from a graph-theoretical perspective. We unified various kinds of specific FHOC problems into two general constrained optimization problems, which can incorporate time-variant costs and a diverse range of constraints. Then, as a central contribution of this study, we established equivalence between general FHOC problems and the SP problem in particular graphs. Two efficient algorithms were afterwards designed to find such an SP. As shown by both time complexity analysis and numerical experiments, our approach can handle the most general problems while maintaining a competitive advantage in computational efficiency. Finally, we note that all SP problems in Problem 1 and 2 can be technically solved by a single SP algorithm, like Dijkstra’s algorithm, though we proposed two custom algorithms for efficiency purpose. That’s why we consider our graph-theoretical approach as a unified framework, which is characterized by high computational efficiency and methodological consistency across a wide range of FHOC problems. Due to the discrete and deterministic nature of BCNs, we believe it is a promising research direction to hybridize the newly developed ASSR with the classical graph theory for more studies on BCNs beyond FHOC. One future work is to adapt this graph-theoretical approach to infinite-horizon optimal control problems.