张安国，李露露, 黄 迟，卢剑权

(1.合肥工业大学 数学学院， 安徽 合肥 230009；2.西南财经大学 经济信息工程学院， 四川 成都 610074；3.东南大学 数学学院，江苏 南京 210096)

0 引言

近年来，随着系统生物学的发展，对基因调控网络的研究吸引了广泛的关注。到目前为止，相关领域的学者们已经提出了大量的数学模型来研究基因调控网络，例如常微分方程，贝叶斯网络，布尔网络等等。布尔网络是由Kauffman教授于1969年提出来的[1]，它对于描述、分析和模拟基因调控网络具有重要作用。布尔网络是由一些离散时间变量构成的，并且每个变量只有两种状态(真和假)，每个变量通过一个逻辑函数更新自己的状态。由于布尔网络结构简单，对于布尔网络的研究在近几年受到了广大学者的密切关注。最近，程代展教授提出了一种新的数学方法——矩阵半张量积，该方法对于布尔网络的研究具有重要促进作用[2-5]。为了进一步研究布尔网络的性质，网络中考虑了控制输入，于是产生了布尔控制网络。基于矩阵半张量积，许多有关布尔网络的问题都已经得到了研究，例如镇定性[6，7]，干扰解耦[8，9]，优化问题[10]，同步[11，12]以及输出跟踪问题[13]等等。

延时在现实中很常见。 例如，人们发现几种外部环境条件可导致细菌生长的时间延迟。 因此，在布尔网络模型中考虑时滞是合理的。对于时滞布尔控制网络的研究已经得到了许多结果。文献[14]和[15]分别研究了时滞布尔控制网络的镇定性和能控性问题；在文献[16]中，作者通过设计状态反馈控制研究了时滞布尔控制网络的同步问题；文献[17]研究了时滞布尔控制网络的输出跟踪问题，并且给出了输出跟踪问题可解的充分必要条件；文献[18]通过设计牵引控制实现了时滞布尔控制网络的稳定性；文献[19]讨论了带有随机时滞的布尔控制网络的稳定性和镇定性问题。

在大量实际系统中，需要根据控制输入将系统的输出驱动到某个理想状态。由于测量技术的限制，不可能直接测量所有状态。然而，可测量的输出可以反映状态。于是，研究布尔网络的输出跟踪问题是有重要现实意义的。文献[20]研究了布尔网络的输出跟踪问题；文献[21]讨论了带有脉冲效应的布尔网络的输出跟踪问题；在文献[13]中，通过设计状态反馈控制，作者讨论了布尔网络的输出追踪给定参考输出信号的问题；在文献[22]中，陈秉权等作者进一步研究了概率布尔控制网络的渐近输出跟踪问题。

关于输出跟踪问题的求解，最常用的控制为状态反馈控制[17]。但是采用状态反馈控制器在求解输出跟踪问题时，一个重要的弊端就是系统在实现输出跟踪之后控制序列可能还会随着状态的更新而发生改变，这将导致极大的资源消耗。一个自然的问题是：能否设计一种控制器使得系统在实现输出跟踪之后控制序列保持不变呢？考虑到时滞在现实生活中的广泛存在性，本文通过设计输出反馈控制器来研究时滞布尔控制网络的输出跟踪问题。

本文的难点主要包含两点：(1) 由于本文考虑的是使用输出反馈求解输出跟踪问题，如何将输出跟踪问题转化为集合镇定问题是本文的第一个难点。(2) 在求解输出反馈增益矩阵时，控制序列需要满足什么要求？这个问题是本文的第二个难点。

本文的主要贡献包括两点：(1) 本文通过设计输出反馈控制器研究了时滞布尔控制网络的输出跟踪问题，并且给出了时滞布尔控制网络输出跟踪问题可以通过输出反馈控制求解的充分必要条件。(2) 针对第二个难点，本文提供了一个有效算法求解输出反馈增益矩阵。

1 预备知识

1.1 常用记号

(2)N+：正整数集； ∅：空集。

(3)coli(M) (rowi(M)): 矩阵M的第i列(行)。

(4)In：n维单位矩阵。

(6) 一个n×p维的矩阵M，如果它的每一列都属于集合Δn，则称其为一个逻辑矩阵。所有的n×p维逻辑矩阵构成的集合记作Ln×p。

(8) |M|：集合M的基数。

(10) [a,b]：a与b之间的所有的整数。

1.2 矩阵半张量积的定义和性质

定义1[2]矩阵M∈Lp×q，N∈Lr×s，l为q和r的最小公倍数，M和N之间的半张量积定义为

关于布尔网络的一些性质在[2-4]中有详细阐述了，因此在本文中我们不再一一列举。

下面介绍一个引理，给出如何将布尔网络逻辑形式转换成代数空间表示形式。

引理1[3]假设f(x1,x2,…,xn)是一个关于变量x1,x2,…,xn的逻辑函数，则可以找到唯一的矩阵Mf∈L2×2n，称为函数f的结构矩阵，使得

f(x1,x2,…,xn)=Mfx1x2…xn。

1.3 模型描述

本文考虑如下系统

(1)

其中x(t)=(x1(t),x2(t),…,xn(t))∈Dn，u(t)=(u1(t),u2(t),…,um(t))∈Dm，y(t)=(y1(t),y2(t),…,yp(t))∈Dp分别表示系统(1)在t时刻的状态、控制输入和输出；fi:Dm+n(τ+1)→D,i=1,2,…n和hj:Dn→D,j=1,2,…p表示系统(1)的逻辑函数。

(2)

其中F∈L2n×2m+n(τ+1)和H∈L2p×2n称为系统(2)的结构矩阵。为了方便，本文用x(t,x(0),…,x(-τ))表示初始状态为x(0),…,x(-τ)时系统在t时刻的状态；同理，y(t,x(0),…,x(-τ))表示初始状态为x(0),…,x(-τ)时系统在t时刻的输出。

由于本文考虑的是带有时滞的布尔网络的输出跟踪问题，下面介绍有关输出跟踪问题的几个定义。

定义2 考虑系统(1)，假设其参考输出信号为y*，称系统(1)的输出跟踪问题是可解的，如果对于任意的初始状态x(0),x(-1),…,x(-τ)，都可找到一个控制序列u(0),u(1),u(2),…和一个正整数k使得

y(t,x(0),…,x(-τ))=y*, ∀t≥k

(3)

定义3 集合S*⊆Δ2n(τ+1)称为系统(1)(或(2))的控制不变子集，如果对于任意的状态x(t)x(t-1)…x(t-τ)∈S*，都存在一个控制u(t)使得x(t+1)x(t)…x(t-τ+1)∈S*。

本文将考虑用输出反馈控制求解带有时滞的布尔控制网络的输出跟踪问题。考虑到时滞对系统的影响，假设设计的控制器的逻辑形式为

(4)

其中gi:Dp(τ+1)→D,i=1,2,…,m为逻辑函数。

u(t)=Gy(t)y(t-1)…y(t-τ)，

(5)

其中G∈L2m×2p(τ+1)称为输出反馈增益矩阵。

2 主要内容

(6)

显然，集合Π(α)包含了所有可以产生参考输出信号的状态。不妨假设集合Π(α)={w1,w2,…,wc}。

下面考虑y(t)y(t-1)…y(t-τ)与x(t)x(t-1)…x(t-τ)之间的关系。通过代数表达式y(t)=Hx(t)，可以得到

(7)

于是，系统(2)可以改写为

(8)

根据以上讨论，可以用如下命题1来描述系统(2)与系统(8)之间的关系。

根据式(5)，控制器u(t)=Gy(t)y(t-1)…y(t-τ)可以改写为

(9)

根据式(1)和式(9)，我们有

另一方面，令z(t)=x(t)x(t-1)…x(t-τ)，那么z(t+1)与z(t)之间的关系可表示为

这里Q=F(I2nτ+m⊗W[2nτ,2n])(I2m⊗Φ2nτ)∈L2n(τ+1)×2m+n(τ+1)。

通过算法1得到了系统的最大控制不变子集S*，并且集合S*中的每个状态都可以产生参考输出信号。下面将根据集合S*给出系统输出跟踪问题在输出反馈控制器下可解的充分必要条件。

(10)

算法1 根据参考输出信号求解系统的最大控制不变子集输入 集合Θ和控制u*=δvβ2p,Sk=∅,k=0,S0=Θ输出 S*While |Sk|>0For any δω12n…δωτ+12n∈Sk If Q-δω12n…δωτ+12n∉Θ Then Sk+1=Θ{δω12n…δωτ+12n} k=k+1 End if End for If Sk+1=SkReturn S*=SkEnd while

关于交运算(∩)，其定义如下,假设X=(xij)，Y=(yij)，定义X∩Y=(xij∩yij)，其中xij∩yij=min{xij,yij}，min{xij,yij}表示元素xij和yij的最小值。

(11)

假设参考输出跟踪问题可解，对于给定的初始状态z(0)，都存在一个状态反馈控制u(z(0))，使得z(0)在控制u(z(0))作用下能够被牵引到集合S*中。因为满足这一要求的控制可能不是唯一的，用U(z(0))表示所有能够将z(0)牵引到集合S*中的状态反馈控制构成的集合，即u(z(0))∈U(z(0))。

关于U(z(0))的求解过程如下。

首先我们通过集合S*构造集合R1(S*)，R2(S*)，…，

R1(S*)={a:Qua∈S*},Rk(S*)={a:Qua∈Rk-1(S*)},∀k≥2。

如果输出跟踪问题可解，那么对于任意的状态a都可以找到唯一的一个k*(k*≥0)使得a∈Rk*(S*)，即存在一个控制u使得Qua=b∈Rk*-1(S*)(定义R0(S*)=S*)。因为满足条件的状态b可能不是唯一的，对于每一个满足条件的状态b，都有一个控制输入u与之对应。于是我们构造集合

U(a)={u:Qua∈Rk*-1(S*)}。

在这里，U(a)表示所有的能够将状态a牵引到集合S*的控制输入构成的集合。

(12)

下面通过定理2给出参考输出跟踪问题能够用输出反馈控制求解的充分必要条件。

定理2 考虑系统(1)，参考输出跟踪问题可以通过输出反馈控制(5)求解当且仅当

uφ≠∅, ∀φ∈{1,…,2p(τ+1)}。

证明(充分性) 对于任意的φ，考虑Π(φ)。当 Π(φ)=∅时，uφ≠∅显然成立。下面考虑Π(φ)≠∅的情形。假设Π(φ)={i1,i2,…,ik}。因为uφ≠∅，假设uφ=u*。通过uφ的构造可知集合Π(φ)中的所有状态都可以在控制u*的作用下牵引到集合S*中。由于输出φ的任意性，于是系统的参考输出跟踪问题可解。

综合以上讨论，我们可以设计如下的算法2来求解输出反馈增益矩阵。

算法2

(1) 首先计算Π(φ)(φ∈[1,2p(τ+1)])；

(2) 通过集合S*，以及代数表达式(2)计算U(z(0)) (z(0)∈Δ2n(τ+1))；

(4) 当Π(φ)≠∅，判断uφ是否为空集。若为空集，则输出跟踪问题不能用输出反馈控制求解。否则输出反馈控制增益矩阵可设计为

G=δ2m[u1，u2，…，u2p(τ+1)]。

注1 通过分析(12)与定理2可知，如果系统可以用输出反馈进行求解，那么它一定也可以通过状态反馈求解。但是如果系统可以通过状态反馈进行求解，它不一定能够通过输出反馈进行求解。

注2 本文主要通过设计输出反馈控制器研究了时滞布尔网络的输出跟踪问题。在文献[13]和[17]中，作者都考虑了系统追踪固定参考输出信号的问题。与文献[13]和[17]相比，本文的主要贡献在于设计了输出反馈控制器而不是状态反馈控制器。另外，因为考虑的都是系统追踪固定参考输出信号的问题，从(5)中可以看出，本文在使用输出反馈控制器来研究系统的输出跟踪问题时，控制输入在一定时间之后将保持不变，这将极大的减少了由于切换控制输入而导致的资源消耗。文献[20]也是通过控制不变子集求解输出跟踪问题。与文献[20]相比，在本文中，控制不变子集中每个状态所对应的控制输入都是相同的。

注3 由于本文考虑的是带有时滞的布尔网络的输出跟踪问题，这将导致在计算的过程中计算复杂度会增加。在求解控制不变子集的时候，时滞的影响将会导致求解过程更加困难。在算法1中，我们通过集合Θ求解了系统的控制不变子集。用a来表示集合Π(y*)中元素的个数，那么集合Θ中元素的个数为aτ+1。另一方面，当u*固定时，第三行执行的次数为O(aτ+1)。事实上，u*是未知的，并且u*可能的个数为2m，所以算法1 的复杂度不超过O(2m×aτ+1)。对于算法2，在第一步中，因为输出状态有2p(τ+1)个，所以第一步的复杂度为O(2p(τ+1))。在第二步中，因为有2n(τ+1)个状态，需要对每个状态分别求可能的控制输入，所以第二步总的时间复杂度为O(2n(τ+1)×2m)。第三步和第四步中，对每个输出状态求解反馈控制，其时间复杂度为O(2p(τ+1))。因此，算法二的复杂度为O(2×2p(τ+1)+2n(τ+1)×2m)。

注4 通过对本文主要定理和算法的分析可以看出：节点的数量，时滞的大小，输出状态的维数和控制输入的维数都会对系统的复杂度产生影响。本文提出的算法一般而言适合研究节点数和时滞都比较小的情况(通常节点数要小于10)。

3 算例分析

考虑如下时滞布尔控制网络的输出跟踪问题

这里x1,x2,y1,y2和u分别表示mRNA, 细胞内乳糖，3号活性胱天蛋白酶(C3a)，8号活性胱天蛋白酶(C8a)，以及α肿瘤坏死因子(TNF)的浓度水平。

设x(t)=x1(t)x2(t)，y(t)=y1(t)y2(t)，z(t)=x(t)x(t-1)，那么上述系统可以改写为

注：横轴表示时间，纵轴表示输出状态中1 所在的位置。图1 初始状态产生输出随时间的变化图

下面计算Rk(S*)。通过计算得到

系统输出随时间的演化情况如图1。

4 总结

本文主要研究了时滞布尔控制网络的输出跟踪问题。首先，利用矩阵半张量积，得到了时滞布尔控制网络的代数形式。其次，通过研究其代数表达式，得到了时滞布尔控制网络输出跟踪问题可解的充要条件。另外，通过研究时滞布尔网络输出跟踪问题可解的充分必要条件，设计了一个算法来求解输出反馈控制器。最后，通过一个例子说明理论结果的可行性。在未来，如何通过设计基于输出反馈的事件触发控制器来实现输出跟踪是一个值得思考的问题。此外，利用矩阵半张量积求解输出跟踪问题，如何降低算法的复杂度也是今后研究的一个重点。