求解大规模线性时不变系统的最优2模型降阶问题的共轭梯度法*

曾泰山，鲁春元，陈 剑

(1.中山大学数学与计算科学学院，广东广州510275;2.广东药学院医药信息工程学院，广东广州510006;3.佛山科学技术学院数学系，广东佛山528000)

曾泰山1，鲁春元2，陈 剑3

(1.中山大学数学与计算科学学院，广东广州510275;2.广东药学院医药信息工程学院，广东广州510006;3.佛山科学技术学院数学系，广东佛山528000)

针对最优2模型降阶问题，提出了适合大规模多输入多输出系统的共轭梯度法。该方法仅需利用一阶导数信息，存储量少，计算复杂度低，且具有超线性收敛性。实验结果显示了算法的有效性。

模型降阶;共轭梯度法;Grassmann流形;线性时+不变系统

近年来，模型降阶越来越受到重视。它可以大大降低大规模系统模拟和控制中的时间复杂度，在超大规模集成电路设计，实时控制，天气预报等领域有着重要应用。关于模型降阶的综述，参见文献[1-2]。

给定矩阵 A∈Rn×n，B∈Rn×p和 C∈Rq×n，线性时不变系统可以表示成:

其中t≥0表示时间变量，u(t)∈Rp表示输入，y(t)∈Rq表示输出，x(t)∈Rn表示系统的状态。n是系统的阶。p和q是系统的输入和输出个数。当系统的阶n很大的时候，由于需要大规模计算和存储，这使得系统的数值模拟和控制非常困难。解决这一问题的一种关键方法是构造一个低阶的系统去逼近原始的高阶系统，并保持相应的系统特性。

由于正交投影能更好地保持系统的稳定性，本文用正交投影的办法构造降阶系统。假设m=n。基于正交投影的模型降阶是选取一个正交矩阵U∈Rn×m，使得UTU=I。我们记A^=UTAU，B^=UTB，C^=CU。通过矩阵U的投影，我们可以获得如下的降阶系统:

定义全阶系统(A，B，C)的传递函数为

方程(3)和(4)描述的降阶系统的传递函数为传递函数G(s)的2范数的平方定义为矩阵积分的迹[9]

定义代价函数

其中，Om表示所有m×m正交矩阵构成的正交群，U∈Rn×m满足UTU=I。由定义可知，等价类[U]中的矩阵的列向量张成的线性子空间相同。在实际计算中，我们将选择其中一个正交矩阵U∈Rn×m来代表整个等价类。关于Grassmann流形的几何性质，参看文献[10-11]。

由方程(6)以及代价函数J(U)的定义 (7)可以看出

也就是说，代价函数J(U)只依赖于U的列向量张成的线性子空间。这意味着代价函数J(U)实际上是Grassmann流形上的光滑函数。因此，最优2模型降阶问题可以改写为Grassmann流形上的最小化问题

通过利用Grassmann流形的几何性质，我们将提出求解问题(9)的数值算法。

在(6)定义的误差传递函数Ge(s)可改写为

实际上，矩阵Ae，Be和Ce定义了一个系统实现为(Ae，Be，Ce)的误差系统。Ge称为是误差系统的传递函数。误差系统的可控性Gramian矩阵Ec和可观测性Gramian矩阵Eo由如下的Lyapunov方程所定义

我们将可控性Gramian矩阵Ec和可观测性Gramian矩阵Eo做如下分划co

Lyapunov方程(10)和(11)可以转化为

代价函数 J(U)可以用 Gramian矩阵 Ec表示如下[9]

我们定义n×m的矩阵JU为代价函数J(U)相对于矩阵U∈Rn×m里面的元素的偏导数:i=1，2，…，n，j=1，2，…，m 。由文献 [7]，偏导数JU可以由以下方式给出:假设P，Q，X和Y是方程(14)、(15)、(16)和 (17)的解。假设R定义为

代价函数J(U)在点U的偏导数为

[11]，在点[U]∈Gr(n，m)处的切空间 T[U]Gr(n，m) 可以表示为

代价函数 J的梯度 ▽J∈ T[U]Gr(n，m) 为

2 共轭梯度法

下面我们给出共轭梯度算法的概要。假设(A，B，C)是全阶系统的状态空间实现。记Uk为第k步获得的模型降阶投影矩阵，k=0，1，…。在共轭梯度算法中，通过以Uk为起点，以Fk为方向的测地线上搜索得到下一个投影矩阵Uk+1。流形上的测地线流形上两点最短路径。以Uk为起点，方向为Fk的Grassmann流形的测地线为

其中 Fk=WkΛkVTk是 Fk奇异值分解，Wk∈ Rn×m，Vk∈ Rm×m，Λk=diag(λ1k，λ2k，…，λmk) 。选取一个合适的步长tk≥0，我们可以构造下一步的模型降阶投影矩阵Uk+1，

为了计算下一步的搜索方向Fk+1，首先计算在点Uk+1处的梯度Gk+1=▽J(Uk+1)。通过求解以下的Sylvester方程来计算矩阵Pk+1，Qk+1，Xk+1和Yk+1，

矩阵Rk+1为

在点Uk+1偏导数JUk+1为JUk+1=2Rk+1。代价函数J在点[Uk+1]∈Gr(n，m)的梯度为

下面将介绍共轭梯度法的搜索方向构造。计算旧搜索方向Fk的平行移动

选取新的搜索方向为共轭梯度方向，即旧搜索方向的平行移动和新的梯度的线性组合

在式子 (28)中，τFk为切向量Fk的沿着测地线的平行移动。切向量在流形上的平行移动是欧氏空间中向量平移的推广。类似于欧氏空间中共轭梯度法，组合系数γk可以通过下面的式子 (Polak-Ribiere法则)得到

其中τGk是梯度Gk的平行移动

值得指出的是，如果令γk=0，则Fk+1=－Gk+1，即搜索方向为负梯度方向，此时共轭梯度法就退化为文献[7]中的梯度流算法。

下面是本文提出的共轭梯度法 (CGA)的主要框架。

Algorithm 1:共轭梯度法 (CGA)

1) 初始化:选取U0∈Rn×m，UT0U0=I，计算G0=▽J(U0)。置F0=－G0;

2)For k=0，1，2，…，N － 1

a)计算Fk的奇异值分解Fk=WkΛkVTk;

b)求解极小化问题

其中:Uk(t) =UkVkcos(tΛk)VTk+Wksin(tΛk)VTk;

令 tk=tmin，Uk+1=Uk(tk);

c)计算梯度Gk+1=▽J(Uk+1);

d)平行移动切向量Fk和Gk到点[Uk+1]:

计算新的搜索方向

3) 获得投影矩阵U=UN，计算

在上述算法中，需要求解极小化问 (30)，这可以通过不完全线性搜索算法来求得，例如Armijo搜索方法，可参见文献[10]。

对于Grassmann流形上的最优化问题，共轭梯度法在满足一定的条件下达到超线性收敛[10-11]。因此，本文提出的共轭梯度算法也能达到超线性收敛。

另外，由于本文所提出的共轭梯度算法不需要求解大规模的 Lyapunov方程，因此计算量大大减少。下面，我们来分析共轭梯度算法的计算复杂性。在本文中，将以乘法次数来衡量计算复杂性。记Ns极小化问题 (30)的最大搜索步数。记N为最大迭代次数。

定理1假设A∈Rn×n是个具有N(A)个非零元素的稀疏矩阵。假设存在一个只需要O(rn+N(A))次乘法运算的线性方程求解方法求解方程其中r是个固定的整数。则共轭梯度法CGA计算复杂度为证明 在第k步迭代中，k=0，1，…，算法CGA的计算花费主要来自以下三个方面:

第一部分是Pk、Qk、Xk和Yk的计算，需要求解Sylvester方程 (23)-(26)。在文献 [7]中指出，如果存在一个只需要O(rn+N(A))次乘法运算的方法求解方程(31)，则利用文献 [12]中的方法求解Sylvester方程 (23)-(26)的计算复杂度只需第二部分是搜索方向Fk的计算。对给定的A，B，C和Uk，需要O(mN(A)+nm2+nmp+nmq)次乘法运算来计算Rk。因此，它需要O(mN(A)+nm2+nmp+nmq)次乘法运算来计算Fk。

第三部分是非精确搜索方法求解最小化问题(30)以获得步长tmin。对极小化问题(30)的每一个搜索步，我们需要计算测地线方程(21)。易知，测地线的计算复杂度为O(nm2)。因为最大的搜索步数为Ns，所以非精确搜索方法求解最小化问题(30)的计算复杂度为O(Nsnm2)。

总而言之，算法CGA每一步迭代需要O(nmr+mN(A)+Nsnm2+nmp+nmq)次乘法次数。由于算法CGA的最大迭代次数为N，因此算法的总的复杂度为如果降阶系统的阶m远远小于原始大规模系统的阶n，且系统的输入输出个数p和q也远远小于n，此时共轭梯度法的计算复杂度为线性O(Nn)，其中N为迭代次数。

3 数值算例

在本节，我们测试比较了本文提出的共轭梯度法的有效性。

例1考虑在区域Ω =(0，1)2上的热传导方程。热传导方程可以有如下形式

其中 u=u(t，x，y) ，(x，y) ∈ Ω ，t∈[0，∞) 。假设微分方程在空间域上等距剖分，格点数为d×d。导出的刚度矩阵A∈Rd2×d2是稀疏的、稳定的，带宽为d。系统的阶为n=d2。假设b1∈Rn是一个所有元素都为1的向量。b2∈Rn是个随机向量。假设B=[b1，b2]，C=BT。此时，构造的系统(A，B，C)是个多输入多输出系统。

我们将使用本文提出的共轭梯度法 CGA将此大规模多输入多输出系统降阶。我们将与以下的方法比较我们的结果:平衡截断方法 BTM[9];快速梯度流算法 FGFA[7]。表1给出了数值结果的比较，2相对误差定义为表中(n，m)给出了系统的阶，其中n是原始大规模系统的阶，m是降阶系统的阶。从表中可以看到，共轭梯度法的计算结果最好。

表1 2相对误差比较Table 1 Comparison of Relative 2error

表1 2相对误差比较Table 1 Comparison of Relative 2error

BTM FGFA CGA(n，m)6.08 e-3 4.05 e-3 4.03 e-3(1600，3) 1.22 e-2 5.96 e-3 5.88 e-3(3600，3)(900，3)1.17 e-2 7.25 e-3 7.10 e-3

为了观察算法的收敛情况，图1给出n=900的时算法FGFA和CGA的收敛曲线。在每一步迭代中，我们计算了2相对误差。从图1中可以看到，两个算法FGFA，CGA的相对误差都随着迭代下降，而算法CGA的收敛速度达到了超线性。

图1 收敛曲线比较Fig.1 Comparison of convergence curves

4 总 结

参考文献:

[1]ANTOULAS A C．Approximation of large-scale dynamical systems[M]．Adv Des Control 6 SIAM，Philadelphia，2005．

[2]GUGERCIN S，ANTOULAS A C．A survey of model reduction by balanced truncation and some new results[J]．International Journal of Control，2004，77(8):748–766．

[3]HUANG X X，YAN W Y，TEO K．Linear-optimal model reduction [J]．IEEE Trans Automat Contr，2001，46:1279–1284．

[4]HYLAND D C，BERNSTEIN D S．The optimal projection equations for model reduction and the relationships among the methods of wilson，skelton and moore[J]．IEEE Trans Automat Contr，1985，30:1201–1211．

[5]LEPSCHY A，MIAN G A，PINATO G，et al．Rational approximation:a non-gradient algorithm[C]//Proc 30th IEEE CDC，1991:2321–2323．

[6]YAN W Y，LAM J．An approximate approach to optimal model reduction[J]．IEEE Trans Automat Contr，1999，44:1341–1358．

[8]GUGERCIN S，ANTOULAS A．C，BEATTIE C．2model reduction for large-scale linear dynamical systems[J]．SIAM J Matrix Anal Appl，2008，30:609–638．

[9]ZHOU K，DOYLE J C，GLOVER K．Robust and optimal control[M]．New Jersey:Prentice-Hall，1996．

[10]ABSIL P A，MAHONY R，SEPULCHRE R．Optimization algorithms on matrix manifolds[M]．Princeton U-niversity Press，2008．

[11]EDELMAN A，ARIAS T A，SMITHS T．The geometry of algorithms with orthogonality constraints[J]．SIAM J Matrix Anal Appl，1999，20:303–353．

[12]VASILYEV D，WHITE J．A more reliable reduction algorithm for behavior model extraction[C]//Proc Int Conf on Computer，Aided Design(ICCAD)，2005:812–819．

A Conjugated Gradient Algorithm for Optimal2Model Reduction of Large Scale Dynamical Systems

ZENG Taishan1，LU Chunyuan2，CHEN Jian3
(1.School of Mathematics and Computational Science，Sun Yat-sen University，Guangzhou 510275，China;2.College of Medical Information Engineering，Guangdong Pharmaceutical University，Guangzhou 510006，China;3.Department of Mathematics，Foshan University，Foshan 528000，China)

A conjugated gradient algorithm with super-linear convergence which is suitable for the optimal2model reduction of the multi-input multi-output large scale dynamical systems is proposed．The proposed algorithm computes only first-order derivative of the cost function．It has low storage requirement and computational cost．Numerical example demonstrates the approximation accuracy and computational efficiency．

model reduction;conjugated gradient method;Grassmann manifold;linear time invariant system

O241.8

A

0529－6579(2011)02－0001－05

2010－04－30

国家自然科学基金资助项目 (10771224);中山大学985项目专项基金资助项目

曾泰山 (1981年生)，男，博士生;E-mail:ztszsu@163.com