一类无穷时域上布尔网络最优控制问题的求解方法

A Method for Solving a Class of Optimal Control Problem of Boolean Networks on Infinite Horizon

摘要: 布尔网络是描述基因调控、细胞分化等系统生物学的有力工具。布尔网络系统的最优控制问题已经成为当前控制领域的研究热点问题之一。本文主要研究布尔网络系统取目标泛函最大化的最优控制问题。首先，给出有限时域上问题的求解方法；其次，在该方法的基础上，对无穷时域上的问题进行研究，给出无穷时域上目标泛函最大化的布尔网络最优控制问题的可解性和求解方法；最后，将方法应用到一个具体的实例上。

Abstract: Boolean network is a powerful tool for describing biological system, such as gene regulation and cell differentiation. The optimal control problem of Boolean networks system has become one of the hot research issues in the field of control on the current. In this paper, we mainly study the optimal control problem of the Boolean network system to maximize the target functional. Firstly, we give the solution method of the problem in finite horizon; secondly, based on the method, the problem in infinite horizon is studied, as well as we give the solvability and solution method of the problem in infinite horizon; in the end, we apply our method to a concrete example.

1. 引言

上世纪六十年代，Jacob和Monod两位教授发现在任何细胞中都包含着几个类似开关一样的调节基因，它们能够打开或者关闭其他的基因 [1] 。1969年，美国学者Kauffman在文献 [1] 的研究基础上，首次提出用布尔网络来刻画细胞和基因调控网络的理论，将基因的表达与不表达用“1”和“0”来表示 [2] 。随后，布尔网络被广泛应用于描述人工智能、神经网络、基因调控等系统，成为系统生物学家、物理学家和系统科学家们共同关心的热点问题。

动态布尔网络系统各节点状态的变化可描述为布尔状态变量的离散时间系统，每一个节点在下一时刻的状态由它及其相邻节点的当前状态确定，数学上可以用布尔函数来描述。为了使用控制论和最优控制理论解决布尔控制系统问题，其关键是如何将逻辑运算表示为相应的代数运算，矩阵半张量积的引入较好地解决了这个瓶颈问题。

矩阵半张量积是中国科学院程代展教授在1998年首次提出的，并在此后多年逐步发展和不断完善。早期的结果被收录在文献 [3] [4] 中，主要结果出版在专著 [5] [6] 中。矩阵半张量积的方法可以将复杂的逻辑推理过程转变为简单的代数形式，该方法激发了国内外众多学者在布尔网络方面的大量研究工作。

关于布尔网络的最优控制问题，2011年，Laschov和Margaliot教授基于矩阵半张量积方法，研究单输入和多输入布尔控制网络Mayer型最优控制问题 [7] [8] ，导出了这一最优控制问题的Pontryagin极大值原理。同年，赵寅等在文献 [9] 中研究了无穷时域上的动态布尔网络控制系统，以及如何寻找最大化序列使其在平均支付条件下目标泛函最大化的问题，运用不动点及极限环等拓扑性质，证明了最优控制可以在乘积空间的一个环上得到，并且具备周期性。李海涛等在文献 [10] 中运用输入–状态关联矩阵研究布尔控制网络Mayer型最优控制问题的求解方法，给出了该问题的最优控制策略。2014年，Fornasini和Elena教授提出求解目标泛函取最小值的逻辑动态网络系统最优控制的一种新的求解方法。在此基础上研究了无穷时域上最优控制问题可解的充分必要条件，并说明可选择一种输入，使其状态经过有限步后，状态–控制对呈现周期性和目标泛函效用为零，因而将无穷时域的最小控制问题转换为有限时域的最优控制问题，从而导出了无穷时域目标泛函最小化问题求解的一种新方法 [11] 。但该方法不能直接用于布尔控制网络目标泛函极大化问题的求解，因其是居于目标泛函非负而设计的算法，不能直接套用经典的最优控制问题中目标泛函的最大化和最小化通过乘以−1相互转化的办法得到。因此，本文主要研究布尔控制网络目标泛函最大化的最优控制问题。本文在给出有限时域上布尔网络控制系统目标泛函最大化求解方法的基础上，给出无穷时域上布尔网络控制最大化问题的可解性和求解方法，最后给出一个具体的算例。

2. 预备知识

2.1. 符号说明

· ${ℝ}^n$ 是实数域上的n维向量空间。

· $M_{m\times n}$ 是 $m\times n$ 维实矩阵集合。

· $I_n$ 是n维单位矩阵。

· ${\Delta }_n=\left\{{\delta }_n^i|i=1,2,\cdots ,n\right\}$ ，其中 ${\delta }_n^i$ 为单位矩阵 $I_n$ 的第i列，特别地 ${\Delta }_2=\Delta$ 。

· 矩阵 $L=\left[{\delta }_n^{i_1},{\delta }_n^{i_2},\cdots ,{\delta }_n^{i_m}\right]\in L_{n\times m}$ ，记为 $L={\delta }_n\left[i_1,i_2,\cdots ,i_m\right]$ 。

· 逻辑变量 $x\in \mathcal{D}=\{0,1\}$ 。

· $\otimes$ 是矩阵的Kronecker积。

2.2. 矩阵半张量积的定义及性质

定义2.2.1 [5] ：给定矩阵 $A\in M_{m\times n},B\in M_{p\times q}$ 。A与B的半张量积记成

$A⋉B=\left(A\otimes I_{t/n}\right)\left(B\otimes I_{t/p}\right)$

其中 $t=lcm\left(n,p\right)$ 是 $\left\{n,p\right\}$ 的最小公倍数。

命题2.2.1 [5] ：矩阵半张量积满足

1) 分配律：

$\left(A+B\right)⋉C=A⋉C+B⋉C$

$C⋉\left(A+B\right)=C⋉A+C⋉B$

2) 结合律

$\left(A⋉B\right)⋉C=A⋉\left(B⋉C\right)$

3) 转置

${\left(A⋉B\right)}^{\text{T}}=B^{\text{T}}⋉A^{\text{T}}$

为进一步实现交换运算，定义如下换位矩阵。

定义2.2.2 [5] ：换位矩阵 $W_{\left[m,n\right]}$ 是一个 $mn\times mn$ 维矩阵，它的行由双指标 $\left(I,J\right)$ 依次标注并按索引 $Id\left(I,J;m,n\right)$ 排列，列由双指标 $\left(i,j\right)$ 依次标注并按索引 $Id\left(i,j;n,m\right)$ 排列，位于 $\left(I,J\right)$ 行 $\left(i,j\right)$ 列的元素为

$W_{\left(I,J\right),\left(i,j\right)}=\{\begin{array}{l}1,I=i且J=j\\ 0,其他\end{array}$

特别地，当 $m=n$ 时，记 $W_{\left[n\right]}=W_{\left[n,n\right]}$ 。

通过换位矩阵，可以实现矩阵的伪交换性质。

命题2.2.2 [5] ：给定矩阵 $A\in M_{m\times n}$ 。

1) 设 $X\in {ℝ}^n$ ， $Y\in {ℝ}^m$ 为列向量，则

$W_{\left[n,m\right]}⋉X⋉Y=Y⋉X,W_{\left[m,n\right]}⋉Y⋉X=X⋉Y$

2) 设 $X\in {ℝ}^n$ 为列向量, $Y\in {ℝ}^m$ 为行向量，则

$X⋉Y=Y⋉W_{\left[n,m\right]}⋉X$

3) 设 $Z\in {ℝ}^t$ 为列向量，则

$Z⋉A=W_{\left[m,t\right]}⋉A⋉W_{\left[t,n\right]}⋉Z=\left(I_t\otimes A\right)⋉Z$

4) 设 $Z\in {ℝ}^t$ 为行向量，则

$A⋉Z=Z⋉W_{\left[m,t\right]}⋉A⋉W_{\left[t,n\right]}=Z⋉\left(I_t\otimes A\right)$

命题2.2.3 [6] ：

哑矩阵 $E_d={\delta }_2\left[1212\right]$ ，具有如下性质：对于任意矩阵 $X,Y\in {\Delta }_2$

$E_{d}XY=Y,E_d{W}_{\left[2\right]}XY=X$

命题2.2.3 [6] ：

降阶矩阵 $M_r=\left[\begin{array}{l}\begin{array}{cc}{\delta }_2^1& 0_2\end{array}\\ \begin{array}{cc}{0}_2& {\delta }_2^2\end{array}\end{array}\right]$ ，具有如下性质：对于任意 $X\in {\Delta }_2$ ，

$X^2=M_{r}X$

命题2.2.6 [6] ：

常用逻辑算子：非( $\neg$ )，合取( $\wedge$ )，析取( $\vee$ )，等值( $\leftrightarrow$ )，异或( $\bar{\vee }$ )的结构矩阵为

$\begin{array}{l}{M}_n={\delta }_2\left[2,1\right],M_c={\delta }_2\left[1,2,2,2\right],M_d={\delta }_2\left[1,1,1,2\right],\\ M_e={\delta }_2\left[1,2,2,1\right],M_m={\delta }_2\left[2,1,1,2\right]\end{array}$

3. 问题描述

3.1. 问题的提出

考虑布尔网络控制系统

$\{\begin{array}{l}{x}_1\left(k+1\right)=f_1\left(x_1\left(k\right),x_2\left(k\right),\cdots ,x_n\left(k\right),u_1\left(k\right),u_2\left(k\right),\cdots ,u_m\left(k\right)\right),\\ x_2\left(k+1\right)=f_2\left(x_1\left(k\right),x_2\left(k\right),\cdots ,x_n\left(k\right),u_1\left(k\right),u_2\left(k\right),\cdots ,u_m\left(k\right)\right),\\ ⋮\\ x_n\left(k+1\right)=f_n\left(x_1\left(k\right),x_2\left(k\right),\cdots ,x_n\left(k\right),u_1\left(k\right),u_2\left(k\right),\cdots ,u_m\left(k\right)\right),\end{array}$ (1)

其中 $x_i\left(k\right)\in \mathcal{D},i=0,1,2,\cdots ,n$ 为系统的状态变量， $u_i\left(k\right)\in \mathcal{D},i=0,1,2,\cdots ,m$ 为系统的控制变量， $f_i:{\mathcal{D}}^{n\text{+}m}\to \mathcal{D}$ 为布尔函数。设初始状态 $x\left(0\right)$ ，终端状态 $x\left(K\right)$ 给定。

容许控制集： $u_{ad}=\left\{u:N\to {\mathcal{D}}^m,u_i\left(k\right)\in \mathcal{D},k=0,1,2,\cdots ,K-1\right\}$

目标泛函：

$J={\lambda }^K{\Phi }_f\left(x\left(K\right)\right)+{\displaystyle {\sum }_{k=0}^{K-1}{\lambda }^k\Phi \left(u\left(k\right),x\left(k\right)\right)}$ (2)

其中 $0<\lambda <1,K>0$ 给定， ${\Phi }_f(\cdot ):{\mathcal{D}}^n\to \mathcal{D},\Phi \left(\cdot ,\cdot \right):{\mathcal{D}}^{n+m}\to \mathcal{D}$ 上的有界效用函数。

问题(P)：寻找 $u^{\ast }\in u_{ad},s.t,J\left(u^{\ast }\right)\ge J\left(u\right),\forall u\in u_{ad}$ 。

3.2. 问题(P)的半张量表达形式

引理3.2.1 [6] ：设任意n元布尔函数 $f\left(x_1,x_2,\cdots ,x_n\right):{\mathcal{D}}^n\to \mathcal{D}$ ，存在唯一的逻辑矩阵 $M_f\in {\mathcal{L}}_{2\times 2^n}$ 。

使在向量形式下

$f\left(x_1,x_2,\cdots ,x_n\right)=M_f{⋉}_i^n{\hat{x}}_i$ (3)

其中 ${\hat{x}}_i\in \Delta$ ， $M_f$ 称为函数f的结构矩阵，(3)式称为布尔函数f的代数形式。

对于受控系统(1)，令 $\hat{x}\left(k\right)={⋉}_i^n{\hat{x}}_i\left(k\right),\hat{u}\left(k\right)={⋉}_i^m{\hat{u}}_i\left(k\right)$ ，由引理3.2.1，对于每一个布尔函数 $f_i$ ，可以找到与之相对应的结构矩阵 $M_i$ ，使得

${\hat{x}}_i\left(k+1\right)=M_i⋉\hat{u}\left(k\right)⋉\hat{x}\left(k\right),i=1,2,\cdots ,n$ (4)

把(4)中n个方程相乘，应用矩阵半张量的性质，得到系统(1)的代数表达式

$\hat{x}\left(k+1\right)=L⋉\hat{u}\left(k\right)⋉\hat{x}\left(k\right)$ (5)

其中 $\hat{x}\left(k\right)={⋉}_i^n{\hat{x}}_i\left(k\right)\in {\Delta }_{2^n}$ 是系统的状态变量， $\hat{u}\left(k\right)={⋉}_i^n{\hat{u}}_i\left(k\right)\in {\Delta }_{2^m}$ 是系统的控制变量， $L\in {\mathcal{L}}_{2^n\times 2^{n+m}}$ ，称为受控系统(1)的状态转移矩阵。

引理3.2.2 [6] ：设任意n元伪布尔函数 $f\left({\hat{x}}_1,{\hat{x}}_2,\cdots ,{\hat{x}}_n\right):{\mathcal{D}}^n\to ℝ$ ，存在唯一的行向量 $V_f\in {ℝ}^{2^n}$ 。

使在向量形式下

$f\left({\hat{x}}_1,{\hat{x}}_2,\cdots ,{\hat{x}}_n\right)=V_f{⋉}_i^n{\hat{x}}_i$ (6)

其中 ${\hat{x}}_i\in \Delta$ ， $V_f$ 称为伪布尔函数f的结构向量。

对于受控系统(1)的目标泛函(2)，令 $\hat{x}\left(k\right)={⋉}_i^n{\hat{x}}_i\left(k\right),\hat{u}\left(k\right)={⋉}_i^m{\hat{u}}_i\left(k\right)$ ，由引理3.2.2，目标泛函可转化为如下形式

$\begin{array}{c}J\left(\hat{u}\right)={\lambda }^K{\Phi }_f\left(\hat{x}\left(K\right)\right)+{\displaystyle {\sum }_{k=0}^{K-1}{\lambda }^k\Phi \left(\hat{u}\left(k\right),\hat{x}\left(k\right)\right)}\\ ={\lambda }^K{\Phi }_f^{\text{T}}⋉\hat{x}\left(K\right)+{\displaystyle {\sum }_{k=0}^{K-1}{\lambda }^k{\Phi }^{\text{T}}⋉\hat{u}\left(k\right)⋉\hat{x}\left(k\right)}\end{array}$ (7)

其中 $0<\lambda <1,K>0$ 给定， ${\Phi }_f(\cdot ):{\mathcal{D}}^n\to ℝ,\Phi \left(\cdot ,\cdot \right):{\mathcal{D}}^{n+m}\to M_{2^m\times 2^n}$ 的有界效用函数， ${\Phi }_f^{\text{T}},{\Phi }^{\text{T}}$ 可按如下方式选取：

$\begin{array}{l}{\Phi }_{{}^f}^{\text{T}}={\alpha }^{\text{T}}\hat{x}\left(K\right)\\ =\left[{\alpha }^{\text{T}}⋉{\delta }_{2^n}^1{\alpha }^{\text{T}}⋉{\delta }_{2^n}^2\cdots {\alpha }^{\text{T}}⋉{\delta }_{2^n}^{2^n}\right]\\ =\left[{\Phi }_f\left({\delta }_{2^n}^1\right){\Phi }_f\left({\delta }_{2^n}^2\right)\cdots {\Phi }_f\left({\delta }_{2^n}^{2^n}\right)\right]\in ℝ\end{array}$

其中 $\alpha \in {ℝ}^{2^n},\left\{a_i\in \alpha |a_i=i,i=1,\cdots ,2^n\right\}$

$\begin{array}{l}{\Phi }^{\text{T}}=\hat{u}{\left(k\right)}^{\text{T}}⋉A⋉\hat{x}\left(k\right)\\ =\left[\Phi \left({\delta }_{2^m}^1,{\delta }_{2^n}^1\right)\cdots \left({\delta }_{2^m}^1,{\delta }_{2^n}^{2^n}\right)\left|\cdots \right|\Phi \left({\delta }_{2^m}^{2^m},{\delta }_{2^n}^1\right)\cdots \Phi \left({\delta }_{2^m}^{2^m},{\delta }_{2^n}^{2^n}\right)\right]\\ =\left[{\Phi }_1^{\text{T}}\left|\cdots \right|{\Phi }_{2^m}^{\text{T}}\right]\in M_{2^m\times 2^n}\end{array}$

其中 $A\in M_{2^m\times 2^n},\left\{a_{ij}\in A|a_{ij}=i+j,i=1,\cdots ,2^m,j=1,\cdots ,2^n\right\}$ 。

进而可以得到如下半张量形式的布尔网络最优控制问题：

受控系统： $\hat{x}\left(k+1\right)=L⋉\hat{u}\left(k\right)⋉\hat{x}(k)$

初始状态 $\hat{x}\left(0\right)\in \Delta$ ，终端状态 $\hat{x}\left(K\right)\in \Delta$ 给定。

容许控制集： ${\hat{u}}_{ad}=\left\{\hat{u}\left(k\right):={⋉}_i^m{\hat{u}}_i\left(k\right),{\hat{u}}_i\left(k\right)\in \Delta ,k=0,1,2,\cdots ,K-1\right\}$

目标泛函：

$J\left(\hat{u}\right)={\lambda }^K{\Phi }_f^{\text{T}}⋉\hat{x}\left(K\right)+{\displaystyle {\sum }_{k=0}^{K-1}{\lambda }^k{\Phi }^{\text{T}}⋉\hat{u}\left(k\right)⋉\hat{x}\left(k\right)}$ (8)

问题( $\text{<mover accent="true"> <mo> P </mo> <mo> ^ </mo> </mover>}$ )：寻找 ${\hat{u}}^{\ast }\in {\hat{u}}_{ad},s.t.,J\left({\hat{u}}^{\ast }\right)\ge J\left(\hat{u}\right),\forall \hat{u}\in {\hat{u}}_{ad}$ 。

4. 问题求解

4.1. 有限时域上布尔网络最优控制问题的求解方法

我们注意到对任意时刻的控制输入 $\hat{u}\left(k\right)={\delta }_{2^m}^i$ ，令 $L⋉\hat{u}\left(k\right)=L_i\in {\mathcal{L}}_{2^n\times 2^n}$ ，则可把布尔控制网络系统(5)看成一个布尔切换系统 [12] [13] [14] [15] ，自然地，布尔切换系统的第i个子系统是一个布尔系统。从而得到

$\hat{x}\left(k+1\right)=L_i\hat{x}\left(k\right),i\in \left[1,2^m\right]$ (9)

系统(5)的状态转移矩阵L可以表示为

$L=\left[L_1{L}_2\cdots L_{2^m}\right]$ (10)

下面考虑有限时域 $\left(K=N<\infty \right)$ 时布尔网络最优控制问题的求解方法。

因对任意n维非正向量 $\omega \left(k\right),k\in \left[0,K\right]$ 和状态 $\hat{x}\left(k\right),k\in \left[0,K\right]$ ，成立下列恒等式：

$0=\omega {\left(0\right)}^{\text{T}}⋉\hat{x}\left(0\right)+{\displaystyle {\sum }_{k=0}^{K-1}\left[\omega {\left(k+1\right)}^{\text{T}}⋉\hat{x}\left(k+1\right)-\omega {\left(k\right)}^{\text{T}}⋉\hat{x}\left(k\right)\right]}-\omega {\left(K\right)}^{\text{T}}⋉\hat{x}\left(K\right)$ (11)

因此，由(8)式-(11)式可得

$\begin{array}{l}J\left(\hat{u}\right)=-\omega {\left(0\right)}^{\text{T}}⋉\hat{x}\left(0\right)+\left({\lambda }^K{\Phi }_f^{\text{T}}+\omega {\left(K\right)}^{\text{T}}\right)⋉\hat{x}\left(K\right)\\ +{\displaystyle {\sum }_{k=0}^{K-1}{\lambda }^k{\Phi }^{\text{T}}⋉\hat{u}\left(k\right)⋉\hat{x}\left(k\right)}-{\displaystyle {\sum }_{k=0}^{K-1}\left[\omega {\left(k+1\right)}^T⋉\hat{x}\left(k+1\right)-\omega {\left(k\right)}^{\text{T}}\hat{x}\left(k\right)\right]}\end{array}$ (12)

又因为 $\hat{u}\left(k\right)={\delta }_{2^m}^i,i\in \left[1,2^m\right]$ ，故对任意状态 $\hat{x}\left(k\right),k\in \left[0,K\right]$ ，成立

$\omega {\left(k\right)}^{\text{T}}\hat{x}\left(k\right)=\left[\omega {\left(k\right)}^{\text{T}}\omega {\left(k\right)}^{\text{T}}\cdots \omega {\left(k\right)}^{\text{T}}\right]⋉\hat{u}\left(k\right)⋉\hat{x}\left(k\right)$ (13)

将(13)式代入(12)式得

$\begin{array}{l}J\left(\hat{u}\right)=-\omega {\left(0\right)}^{\text{T}}⋉\hat{x}\left(0\right)+\left({\lambda }^K{\Phi }_f^{\text{T}}+\omega {\left(K\right)}^{\text{T}}\right)⋉\hat{x}\left(K\right)\\ +{\displaystyle {\sum }_{k=0}^{K-1}\left({\lambda }^k{\Phi }^{\text{T}}-\omega {\left(k+1\right)}^{\text{T}}L+\left[\omega {\left(k\right)}^{\text{T}}\cdots \omega {\left(k\right)}^{\text{T}}\right]\right)}⋉\hat{u}\left(k\right)⋉\hat{x}\left(k\right)\end{array}$ (14)

令

$\begin{array}{c}\psi {\left(k\right)}^{\text{T}}=\left[{\psi }_1{\left(k\right)}^{\text{T}}{\psi }_2{\left(k\right)}^{\text{T}}\cdots {\psi }_{2^m}{\left(k\right)}^{\text{T}}\right]\\ =\left[{\lambda }^k{\Phi }_1^{\text{T}}{\lambda }^k{\Phi }_2^{\text{T}}\left|\cdots \right|{\lambda }^k{\Phi }_{2^m}^T\right]-\omega {\left(k+1\right)}^{\text{T}}\left[L_1{L}_2\cdots L_{2^m}\right]+\left[\omega {\left(k\right)}^{\text{T}}\omega {\left(k\right)}^{\text{T}}\cdots \omega {\left(k\right)}^{\text{T}}\right]\end{array}$ (15)

即(14)式可等价表示为

$\begin{array}{l}J\left(\hat{u}\right)=-\omega {\left(0\right)}^{\text{T}}⋉\hat{x}\left(0\right)+\left({\lambda }^K{\Phi }_f^{\text{T}}+\omega {\left(K\right)}^{\text{T}}\right)⋉\hat{x}\left(K\right)\\ +{\displaystyle {\sum }_{k=0}^{K-1}{\left[{\psi }_1{\left(k\right)}^{\text{T}}{\psi }_2{\left(k\right)}^{\text{T}}\cdots {\psi }_{2^m}\left(k\right)\right]}^{\text{T}}⋉\hat{u}\left(k\right)⋉\hat{x}\left(k\right)}\end{array}$ (16)

为使目标泛函(16)式达到最大，根据贝尔曼动态规划方法的最优性原理，我们希望每一步的取值达到最大，即取 $\omega \left(K\right)$ 为零向量时使终端达到最大，对 $\forall k\in \left[0,K-1\right],j\in \left[1,2^n\right]$ ，存在 $i\in \left[1,2^m\right]$ 使得向量 ${\left[{\psi }_i\left(k\right)\right]}_j=0$ ，否则 ${\left[{\psi }_i\left(k\right)\right]}_j\le 0$ 。

因此，我们设计如下算法，使目标泛函(16)式达到最大：

Step 1：设向量 $\omega \left(K\right)=0_{2^n}^{\text{T}}$ ；

Step 2：当 $k=K-1,K-2,\cdots ,1,0$ ，向量 $\omega \left(k\right)$ 的第j项为

$-{\left[\omega \left(k\right)\right]}_j=\underset{i\in \left[1,2^m\right]}{\max }\left({\lambda }^k{\left[{\Phi }_i\right]}_j-{\left[\omega {\left(k+1\right)}^{\text{T}}{L}_i\right]}_j\right),\forall j\in \left[1,2^n\right];$

Step 3：使目标泛函(16)式达到最大值的控制输入序列 $\hat{u}\left(k\right),k\in \left[0,K-1\right]$ 可选择如下：

对 $\forall \hat{x}\left(k\right)={\delta }_{2^n}^j$ ，可得 $\hat{u}\left(k\right)={\delta }_{2^m}^{i^{*}\left(j,k\right)}$ 。

其中

$i^{*}\left(j,k\right)=\arg \underset{i\in \left[1,2^m\right]}{\max }\left({\lambda }^k{\left[{\Phi }_i\right]}_j-{\left[\omega {\left(k+1\right)}^{\text{T}}{L}_i\right]}_j\right),j\in \left[1,2^n\right]$

注1：事实上，根据算法的设计，最优控制序列 $\hat{u}\left(k\right)$ 实际上是一个反馈控制，即

$\hat{u}\left(k\right)={\rm H}\left(k\right)\hat{x}(k)$

其中

${\rm H}\left(k\right)=\left[{\delta }_{2^m}^{i^{*}\left(1,k\right)}{\delta }_{2^m}^{i^{*}\left(2,k\right)}\cdots {\delta }_{2^m}^{i^{*}\left(2^n,k\right)}\right]$

注2：根据贝尔曼动态规划方法，令值函数为

$\{\begin{array}{l}{J}_K\left(\hat{x}\left(K\right)\right)=\underset{\hat{u}\left(K-1\right)\in {\Delta }_{2^m}}{\max }{\lambda }^K{\Phi }_f^{\text{T}}⋉{\delta }_{2^n}^j\\ J_{\left[0,K-1\right]}\left(\hat{x}\left(k\right),\hat{u}(\cdot )\right)=J_{\left[0,K-1\right]}\left({\delta }_{2^n}^j,\hat{u}(\cdot )\right)=\underset{\hat{u}(\cdot )\in {\Delta }_{2^m}}{\max }\left\{{\lambda }^K{\Phi }_f^{\text{T}}⋉{\delta }_{2^n}^j+J_{k+1}\left({\delta }_{2^n}^j,\hat{u}(\cdot )\right)\right\}\end{array}$

由本文设计的算法可知

$-{\left[\omega \left(k\right)\right]}_j=\underset{\hat{u}\left(k\right)\in {\Delta }_{2^m}}{\max }J\left({\delta }_{2^n}^j,\hat{u}(\cdot )\right)$

$-{\left[\omega \left(k\right)\right]}_j=\underset{\hat{u}\left(k\right)\in {\Delta }_{2^m}}{\max }\left\{{\lambda }^K{\Phi }_f^{\text{T}}⋉{\delta }_{2^n}^j+J_{k+1}\left({\delta }_{2^n}^j,\hat{u}(\cdot )\right)\right\}$

$-{\left[\omega \left(k\right)\right]}_j=\underset{\hat{u}\left(k\right)\in {\Delta }_{2^m}}{\max }\left\{{\lambda }^K{\Phi }_f^{\text{T}}⋉{\delta }_{2^n}^j\right\}+\underset{\hat{u}(k)\in {\Delta }_{2^m}}{\max }\left\{J_{k+1}\left({\delta }_{2^n}^j,\hat{u}\left(k\right)\right)\right\}$

从中我们得到，本文给出的方法与动态规划算法的求解结果相一致。

综上，我们得到了有限时域上布尔网络最优控制问题的求解方法。

4.2. 无穷时域上布尔网络最优控制的求解方法

下面我们考虑受控系统(5)支配下，无穷时域 $\left(K\to \infty \right)$ 上的布尔网络最优控制问题。设目标泛函为

$J\left(\hat{u}\right)={\displaystyle {\sum }_{k=0}^{\infty }{\lambda }^k\Phi \left(\hat{u}\left(k\right),\hat{x}\left(k\right)\right)}$ (17)

其中 $0<\lambda <1$ ， $\Phi \left(\cdot ,\cdot \right):{\mathcal{D}}^{n+m}\to M_{2^n\times 2^m}$ 上的有界效用函数。

对于无穷时域上的最优控制问题，其最优控制问题的可解性等价于目标泛函中无穷级数的收敛性。如下定理给出无穷时域上布尔网络最优控制问题的求解方法。

定理4.2.1：对布尔控制系统(5)以及目标泛函(17)，总存在最优控制 ${\hat{u}}^{*}\left(k\right)$ ，且当有限时域中终端步K取得充分大时，从有限时域上求得的最优控制序列 $\left\{\hat{u}{\left(1\right)}^{*},\hat{u}{\left(2\right)}^{*},\cdots ,\hat{u}{\left(K\right)}^{*}\right\}$ 与无穷时域上的最优控制序列的前K个值一致。

证明：令 $\mathcal{M}=\underset{\hat{u}\left(k\right),\hat{x}\left(k\right)\in \Delta }{\max }\left|\Phi \left(\hat{u}\left(k\right),\hat{x}\left(k\right)\right)\right|$ ，则由 $\hat{u}\left(k\right),\hat{x}\left(k\right)$ 都是有限维向量， $\Phi \left(\cdot ,\cdot \right):{\mathcal{D}}^{n+m}\to M_{2^m\times 2^n}$

上有界函数，可知 $\mathcal{M}<+\infty$ ，所以无穷级数 ${\displaystyle {\sum }_{k=0}^{\infty }{\lambda }^k\Phi \left(\hat{u}\left(k\right),\hat{x}\left(k\right)\right)}$ 的优级数为 $\mathcal{M}{\displaystyle {\sum }_{k=0}^{\infty }{\lambda }^k}$ 是收敛的。

即无穷级数 ${\displaystyle {\sum }_{k=0}^{\infty }{\lambda }^k\Phi \left(\hat{u}\left(k\right),\hat{x}\left(k\right)\right)}$ 一致收敛。故对任意给定的 $\epsilon >0$ ，存在足够大的p，使得

$\left|{\displaystyle {\sum }_{k=p}^{\infty }{\lambda }^k\Phi \left(\hat{u}\left(k\right),\hat{x}\left(k\right)\right)}\right|\le {\displaystyle {\sum }_{k=p}^{\infty }{\lambda }^p\mathcal{M}}=\frac{{\lambda }^{p+1}}{1-\lambda }\mathcal{M}<\frac{\epsilon }{2}$

不妨设 $K=p$ ，考虑受控系统如下求解有限时域最优控制问题

$J_K\left({\hat{x}}_0,\hat{u}\right)={\displaystyle {\sum }_{k=0}^K{\lambda }^k\Phi (\hat{u}\left(k\right),\hat{x}(k))}$

若 $\left\{\hat{u}{\left(1\right)}^{*},\hat{u}{\left(2\right)}^{*},\cdots ,\hat{u}{\left(K\right)}^{*}\right\}$ 为 $J_K\left({\hat{x}}_0,\hat{u}\right)$ 的最优控制，若 $\left\{{\hat{u}}^{\infty }{\left(1\right)}^{*},{\hat{u}}^{\infty }{\left(2\right)}^{*},\cdots ,{\hat{u}}^{\infty }{\left(K\right)}^{*},\cdots \right\}$ 为 $J_{\infty }\left({\hat{x}}_0,\hat{u}\right)$ 的最优控制，则有

$J_K\left({\hat{x}}_0;\hat{u}{\left(1\right)}^{*},\hat{u}{\left(2\right)}^{*},\cdots ,\hat{u}{\left(K\right)}^{*}\right)-J_{\infty }\left({\hat{x}}_0;{\hat{u}}^{\infty }{\left(1\right)}^{*},{\hat{u}}^{\infty }{\left(2\right)}^{*},\cdots ,{\hat{u}}^{\infty }{\left(K\right)}^{*},\cdots \right)\le \epsilon$

否则

$\begin{array}{l}J\left({\hat{x}}_0;\hat{u}{\left(1\right)}^{*},\hat{u}{\left(2\right)}^{*},\cdots ,\hat{u}{\left(K\right)}^{*},{\hat{u}}^{\infty }{\left(K+1\right)}^{*},\cdots \right)-J\left({\hat{x}}_0;{\hat{u}}^{\infty }{\left(1\right)}^{*},{\hat{u}}^{\infty }{\left(2\right)}^{*},\cdots ,{\hat{u}}^{\infty }{\left(K\right)}^{*},\cdots \right)\\ ={\displaystyle {\sum }_{k=0}^K{\lambda }^k\left[\Phi \left(\hat{u}{\left(k\right)}^{*},\hat{x}{\left(k\right)}^{*}\right)-\Phi \left({\hat{u}}^{\infty }{\left(k\right)}^{*},{\hat{x}}^{\infty }{\left(k\right)}^{*}\right)\right]}\\ +{\displaystyle {\sum }_{k=K}^{\infty }{\lambda }^k\left[\Phi \left(\hat{u}{\left(k\right)}^{*},\hat{x}{\left(k\right)}^{*}\right)-\Phi \left({\hat{u}}^{\infty }{\left(k\right)}^{*},{\hat{x}}^{\infty }{\left(k\right)}^{*}\right)\right]}\\ >\epsilon -2\cdot \frac{\epsilon }{2}=0\end{array}$

这与 ${\hat{u}}^{\infty }{\left(k\right)}^{*}$ 为最优控制相矛盾，即无穷时域上的最优控制存在。 $\square$

注3：当K取得足够大时，无穷时域上的最优控制序列 $\left\{{\hat{u}}^{\infty }{\left(1\right)}^{*},{\hat{u}}^{\infty }{\left(2\right)}^{*},\cdots ,{\hat{u}}^{\infty }{\left(K\right)}^{*},\cdots \right\}$ 可从有限时域得到的最优控制序列 $\left\{\hat{u}{\left(1\right)}^{*},\hat{u}{\left(2\right)}^{*},\cdots ,\hat{u}{\left(K\right)}^{*}\right\}$ 逐次逼近得到。

5. 应用实例

上述方法为我们求解布尔网络最优控制问题提供了一种思路清晰，行之有效的求解方法，下面将给出实例具体验证方法的可行性。

考虑如下布尔网络最优控制问题

$\{\begin{array}{l}{x}_1\left(k+1\right)=\neg \left(x_1\left(k\right)\leftrightarrow x_2\left(k\right)\right),\\ x_2\left(k+1\right)=\neg \left(x_2\left(k\right)\leftrightarrow x_3\left(k\right)\right),\\ x_3\left(k+1\right)=u\left(k\right)\wedge x_1\left(k\right),\end{array}$ (18)

其中 $x_i\left(k\right)\in \mathcal{D},i=1,2,3$ 为系统状态变量， $u\left(k\right)\in \mathcal{D}$ 为系统控制变量， $f_i:{\mathcal{D}}^{\text{3+}1}\to \mathcal{D}$ 为布尔函数。设初值 $x_1\left(0\right)=1,x_2\left(0\right)=0,x_3\left(0\right)=1$ ,终端 $x_1\left(4\right)=0,x_2\left(4\right)=0,x_3\left(4\right)=0$ ,目标泛函为

$J={\left(0.9\right)}^4{\Phi }_f\left(x\left(4\right)\right)+{\displaystyle {\sum }_{k=0}^3{\left(0.9\right)}^k\Phi \left(u\left(k\right),x\left(k\right)\right)}$ (19)

其中上的有界效用函数。

问题(Q)求使目标泛函(19)式达到最大的控制序列 $\left\{u^{*}\left(0\right),u^{*}\left(1\right),u^{*}\left(2\right),u^{*}\left(3\right)\right\}$ 。

解：首先，将布尔网络系统的最优控制问题转化代数形式的最优控制问题；

①将上述布尔控制网络系统(18)转化为代数形式

令 $\hat{x}\left(k\right)={\hat{x}}_1\left(k\right)⋉{\hat{x}}_2\left(k\right)⋉{\hat{x}}_3\left(k\right)\in \Delta$ 为系统状态变量， $\hat{u}\left(k\right)\in \Delta$ 为系统控制变量

$\{\begin{array}{l}{\hat{x}}_1\left(k+1\right)=M_m⋉{\hat{x}}_1\left(k\right)⋉{\hat{x}}_2\left(k\right),\\ {\hat{x}}_2\left(k+1\right)=M_m⋉{\hat{x}}_2\left(k\right)⋉{\hat{x}}_3\left(k\right),\\ {\hat{x}}_3\left(k+1\right)=M_c⋉\hat{u}\left(k\right)⋉{\hat{x}}_1\left(k\right),\end{array}$ (20)

$\begin{array}{l}\hat{x}\left(k+1\right)=M_m⋉{\hat{x}}_1\left(k\right)⋉{\hat{x}}_2\left(k\right)⋉M_m⋉{\hat{x}}_2\left(k\right)⋉{\hat{x}}_3\left(k\right)⋉M_c⋉\hat{u}\left(k\right)⋉{\hat{x}}_1\left(k\right)\\ =M_m⋉\left(I_4\otimes M_m\right)⋉\left(I_2\otimes M_r\right)⋉\left(I_8\otimes M_c\right)⋉W_{\left[4,8\right]}⋉\left(I_2\otimes M_r\right)⋉\hat{u}\left(k\right)⋉\hat{x}\left(k\right)\\ =L⋉\hat{u}\left(k\right)⋉\hat{x}(k)\end{array}$

此时，布尔控制网络系统的状态转移矩阵为

$L={\delta }_8\left[7513426886244268\right]\in {\mathcal{L}}_{2^3\times 2^{3+1}}$

②将目标泛函(19)式中的伪布尔函数转化为半张量形式

$J={\left(0.9\right)}^4{\Phi }_f^{\text{T}}⋉\hat{x}\left(4\right)+{\displaystyle {\sum }_{k=0}^3{\left(0.9\right)}^k{\Phi }^{\text{T}}⋉\hat{u}\left(k\right)⋉\hat{x}\left(k\right)}$ (21)

其中

$\begin{array}{l}{\Phi }_{{}^f}^{\text{T}}={\alpha }^{\text{T}}\hat{x}\left(K\right)\\ =\left[{\alpha }^{\text{T}}⋉{\delta }_8^1{\alpha }^{\text{T}}⋉{\delta }_8^2\cdots {\alpha }^{\text{T}}⋉{\delta }_8^8\right]\\ =\left[12345678\right]\end{array}$

其中 $\alpha \in {ℝ}^8,\left\{a_i\in \alpha |a_i=i,i=1,2,\cdots ,8\right\}$ 。

$\begin{array}{l}{\Phi }^{\text{T}}=\hat{u}{\left(k\right)}^{\text{T}}⋉A⋉\hat{x}\left(k\right)\\ =\left[\Phi \left({\delta }_2^1,{\delta }_8^1\right)\cdots \Phi \left({\delta }_2^1,{\delta }_8^8\right)\Phi \left({\delta }_2^2,{\delta }_8^1\right)\cdots \Phi \left({\delta }_2^2,{\delta }_8^8\right)\right]\\ =\left[{\Phi }_1^{\text{T}}{\Phi }_2^{\text{T}}\right]\\ =\left[\begin{array}{l}23456789\\ 345678910\end{array}\right]\end{array}$

其中 $A\in M_{2\times 8},\left\{a_{ij}\in A|a_{ij}=i+j,i=1,2,j=1,2,\cdots ,8\right\}$ 。

③求解问题( $\text{<mover accent="true"> <mo> Q </mo> <mo> ^ </mo> </mover>}$ )给定初值状态 $\hat{x}\left(0\right)={\delta }_8^3$ ，终端状态 $\hat{x}\left(4\right)={\delta }_8^8$ ，求使受控系统在半张量形式下的目标泛函(21)式达到最大的最优控制序列 $\left\{{\hat{u}}^{*}\left(0\right),{\hat{u}}^{*}\left(1\right),{\hat{u}}^{*}\left(2\right),{\hat{u}}^{*}\left(3\right)\right\}$ 。

其次，求解问题( $\text{<mover accent="true"> <mo> Q </mo> <mo> ^ </mo> </mover>}$ )；

应用4.1.中有限时域上布尔网络最优控制的求解方法。

①将布尔控制系统转化为切换布尔系统

$\hat{x}\left(k+1\right)=L\hat{u}\left(k\right)\hat{x}(k)$

$\hat{x}\left(k+1\right)=L_i\hat{x}\left(k\right),i\in \left[1,2\right]$

$\begin{array}{l}{L}_1=L⋉{\delta }_2^1={\delta }_8\left[75134268\right]\\ L_2=L⋉{\delta }_2^2={\delta }_8\left[86244268\right]\end{array}$

②根据递归算法，可得

$\omega \left(4\right)=0_{2^3}^{\text{T}}=\left[000000000\right]$

$-\omega \left(3\right)=\underset{i\in \left[1,2\right],j\in \left[1,8\right]}{\max }\left({\lambda }^3{\left[{\Phi }_i\right]}_j-{\left[\omega {\left(4\right)}^{\text{T}}{L}_i\right]}_j\right)=\left[2.192.923.654.375.105.836.567.29\right]$

$\Rightarrow {\hat{x}}^{*}\left(3\right)={\delta }_8^8$

${\rm H}\left(3\right)=\left[{\delta }_2^{i^{*}\left(1,3\right)}{\delta }_2^{i^{*}\left(2,3\right)}\cdots {\delta }_2^{i^{*}\left(8,3\right)}\right]=\left[{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2\right]$

$\Rightarrow {\hat{u}}^{*}\left(3\right)={\delta }_2^2$

$-\omega \left(2\right)=\underset{i\in \left[1,2\right],j\in \left[1,8\right]}{\max }\left({\lambda }^2{\left[{\Phi }_i\right]}_j-{\left[\omega {\left(3\right)}^{\text{T}}{L}_i\right]}_j\right)=\left[9.729.076.979.2310.049.4013.1215.39\right]$

$\Rightarrow {\hat{x}}^{*}\left(2\right)={\delta }_8^8$

${\rm H}\left(2\right)=\left[{\delta }_2^{i^{*}\left(1,2\right)}{\delta }_2^{i^{*}\left(2,2\right)}\cdots {\delta }_2^{i^{*}\left(8,2\right)}\right]=\left[{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2\right]$

$\Rightarrow {\hat{u}}^{*}\left(2\right)={\delta }_2^2$

$-\omega \left(1\right)=\underset{i\in \left[1,2\right],j\in \left[1,8\right]}{\max }\left(\lambda {\left[{\Phi }_i\right]}_j-{\left[\omega {\left(2\right)}^{\text{T}}{L}_i\right]}_j\right)=\left[18.0913.0013.5714.6315.5316.2717.5024.39\right]$

$\Rightarrow {\hat{x}}^{*}\left(1\right)={\delta }_8^8$

${\rm H}\left(1\right)=\left[{\delta }_2^{i^{*}\left(1,1\right)}{\delta }_2^{i^{*}\left(2,1\right)}\cdots {\delta }_2^{i^{*}\left(8,1\right)}\right]=\left[{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2\right]$

$\Rightarrow {\hat{u}}^{*}\left(1\right)={\delta }_2^2$

$-\omega \left(0\right)=\underset{i\in \left[1,2\right],j\in \left[1,8\right]}{\max }\left({\left[{\Phi }_i\right]}_j-{\left[\omega {\left(1\right)}^{\text{T}}{L}_i\right]}_j\right)=\left[27.3920.2722.0920.6321.6321.0025.2734.39\right]$

${\hat{x}}^{*}\left(0\right)={\delta }_8^3$ 给定

${\rm H}\left(0\right)=\left[{\delta }_2^{i^{*}\left(1,0\right)}{\delta }_2^{i^{*}\left(2,0\right)}\cdots {\delta }_2^{i^{*}\left(8,0\right)}\right]=\left[{\delta }_2^2{\delta }_2^2{\delta }_2^1{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2{\delta }_2^2\right]$

$\Rightarrow {\hat{u}}^{*}\left(0\right)={\delta }_2^1$

$\begin{array}{c}J\left({\hat{u}}^{*}\right)={\left(0.9\right)}^4{\Phi }_f^{\text{T}}⋉\hat{x}\left(4\right)+{\displaystyle {\sum }_{k=0}^3{\left(0.9\right)}^k{\Phi }^{\text{T}}⋉\hat{u}\left(k\right)⋉\hat{x}\left(k\right)}\\ =0.66\times \left[12345678\right]⋉{\delta }_8^8+\underset{{\hat{u}}^{*}(\cdot )}{\max }J\left({\delta }_8^3,\hat{u}(\cdot )\right)\\ =5.25-\omega \left(0\right)⋉{\delta }_8^3\\ =27.34\end{array}$

最优控制序列为 ${\hat{u}}^{*}\left(0\right)={\delta }_2^1,{\hat{u}}^{*}\left(1\right)={\hat{u}}^{*}\left(2\right)={\hat{u}}^{*}\left(3\right)={\delta }_2^2$ 。

状态轨迹为 ${\hat{x}}^{*}\left(0\right)={\hat{x}}^{*}\left(1\right)={\hat{x}}^{*}\left(2\right)={\hat{x}}^{*}\left(3\right)={\hat{x}}^{*}\left(4\right)={\delta }_8^8$ 。

最后，将问题( $\text{<mover accent="true"> <mo> Q </mo> <mo> ^ </mo> </mover>}$ )代数形式的求解结果返回问题(Q)的逻辑形式，得到

最优控制序列为 $\left\{u^{*}\left(0\right)\to u^{*}\left(1\right)\to u^{*}\left(2\right)\to u^{*}\left(3\right)\right\}=\left\{1\to 0\to 0\to 0\right\}$ 。

状态轨迹为 $\begin{array}{l}\left\{x^{*}\left(0\right)\to x^{*}\left(1\right)\to x^{*}\left(2\right)\to x^{*}\left(3\right)\to x^{*}\left(4\right)\right\}\\ =\left\{{\left(101\right)}^{\text{T}}\to {\left(000\right)}^{\text{T}}\to {\left(000\right)}^{\text{T}}\to {\left(000\right)}^{\text{T}}\to {\left(000\right)}^{\text{T}}\right\}\end{array}$ 。

6. 总结

本文运用矩阵半张量积的方法，将布尔控制系统转化为代数离散时间系统，从而考虑经典的布尔网络最优控制问题，在找到有限时域上使目标泛函达到最大值时的最优控制及状态轨迹方法的基础上，重点研究无穷时域上的最优控制问题，给出无穷时域上布尔网络控制最大化目标泛函问题的可解性和求解方法，并应用到具体的实例，说明该方法可行。

基金项目

本文获得国家自然科学基金项目(11761021)资助。