K-power双线性系统基于Laguerre函数的保结构模型降阶方法.pdf

应用数学MATHEMATICA APPLICATA2023,36(4):929-939K-power双线性系统基于Laguerre函数的保结构模型降阶方法金辉1,肖志华1,祁振中2(1.长江大学信息与数学学院,湖北 荆州 434023;2.西北大学数学学院,陕西 西安 710127)摘要：针对K-power双线性系统,提出一种基于Laguerre函数的保结构模型降阶方法.该方法首先将K-power双线性系统转化为一般的双线性系统,再利用矩阵指数函数的Laguerre函数展开,计算出双线性系统交叉Gram矩阵的近似低秩因子,从而构造K-power双线性系统各子系统相应的投影变换,进而生成该系统的近似平衡系统.然后,通过截断较小的奇异值对应的状态,得到降阶模型.该降阶过程计算高效,且具有一定的自适应性.最后,数值实验验证了算法的有效性.关键词：模型降阶;K-power双线性系统;Laguerre函数;交叉Gram矩阵中图分类号：O231.2AMS(2010)主题分类：41A10;93A15文献标识码：A文章编号：1001-9847(2023)04-0929-111.引言在工程应用领域中,许多物理系统或现象都可以用数学模型来描述,然而大多数的数学模型规模都极其复杂,使得对这些系统的直接仿真模拟和理论分析变得十分困难.因此,有效减小系统规模和缩短计算仿真时间的需求日益上升.模型降阶正是基于此思想,将一个较大规模的复杂系统转化为一个近似的较小系统,同时保留了原始系统的一些重要性能,例如稳定性、结构性和无源性等1.模型降阶方法已经普遍用于超大规模集成电路模拟、计算电磁学、微机电系统和流体力学等领域2.模型降阶自提出以来,已经形成了一系列成熟的降阶理论和方法,例如Krylov子空间方法、平衡截断方法、本征正交分解方法和正交多项式方法等.其中平衡截断是模型降阶最著名的方法之一3,其基本思想是计算出控制系统的可控Gram矩阵和可观Gram矩阵,进一步构造合适的映射子空间获得降阶模型.该方法可以构造出保持原始系统稳定性的降阶模型,并且给出近似先验误差.然而其需求解两个大规模Lyapunov方程,该过程计算量大且计算时间长,这导致了其应用的局限性.在这种情况下,许多学者作出相应改进,提出了近似平衡截断方法4.该方法以一种数值高效的方式获得近似平衡系统,被广泛应用于线性系统及非线性系统的模型降阶中.对于线性系统,相应的模型降阶方法已经日趋完善.而实际的物理系统本质是非线性的,许多物理系统都可以通过双线性系统模型来描述.双线性系统是一类特殊的非线性系统,是连接线性系统和非线性系统之间的桥梁,该系统经常出现在非线性电路、热过程和生态系统等科学领域5.关于双线性系统的模型降阶,已经发展了多种方法,例如平衡相关方法、Krylov子收稿日期：2022-09-21基金项目：国家自然科学基金项目(62273059,61803046)通讯作者：肖志华,男,汉族,湖北人,副教授,研究方向:模型降阶,数值代数.930应用数学2023空间方法和插值法等68.值得注意的是,Shaker和Tahavori在文9中定义了双线性系统的交叉Gram矩阵.该交叉Gram矩阵同时反映了系统的可控性和可观性,并且当双线性系统有界输入有界输出(BIBO)稳定时,其为广义Sylvester方程的解.作为双线性系统的一个特殊子类,K-power双线性系统的应用也很广泛,涉及液压驱动建模、多项式系统、非线性系统辨识等领域1011.关于K-power双线性系统的模型降阶方法,在文10中,Al-Baiyat和Bettayeb首先证明了K-power双线性系统具有块对角可控Gram矩阵和可观Gram矩阵,并提出了一种简化求解Lyapunov方程计算的保结构平衡截断方法.XIE和Syrmos11证明了K-power双系统的双线性广义奇异摄动近似(GSPA)可以转化为-交互系统的直接截断,基于此,提出了一种基于双线性-交互系统分析的方法.WANG和JIANG12还成功地将矩匹配方法和最优H2方法应用于K-power双线性系统的模型降阶,进一步得到保持原始系统结构特征的降阶模型.本文基于Laguerre函数和双线性系统的交叉Gram矩阵,提出了一种K-power双线性系统的保结构模型降阶方法.该方法首先将其等价双线性系统的交叉Gram矩阵用Laguerre函数进行展开,并结合其正交性,进一步给出交叉Gram矩阵的低秩因子,最后通过构造K-power双线性系统各子系统的投影变换得到降阶模型.相比于经典的平衡截断方法,此方法不需要直接求解Sylvester方程来获得交叉Gram矩阵,具有较高的计算效率和灵活度,且具有一定的自适应性.本文首先介绍Laguerre函数的相关重要性质,接着详细给出保结构的K-power双线性系统降阶过程,并给出了相应的降阶算法.最后,通过数值实验有效验证了算法的有效性.2.矩阵指数函数的Laguerre函数展开矩阵指数函数eAt又称状态转移矩阵,在现代控制理论中,可以将其用于状态方程的求解.本节的目标是利用Laguerre函数13对矩阵指数函数进行展开.定义2.1多项式:l0(t)=1,li(t)=eti!didti(tiet),i=1,2,称为Laguerre多项式.定义2.2函数:i(t)=2etli(2t),称为Laguerre函数.其中 0是一个常数,称为Laguerre参数.性质2.1在时间区间0,+),Laguerre函数满足如下正交性14:0i(t)j(t)dt=0,i=j,1,i=j.性质2.2对Laguerre函数作Laplace变换14,可以得到:i(s)=L(i(t)=2s+(s s+)i,i=0,1,.对于一般指数函数eat(Re(a)0),可以得到其Laguerre函数展开式15eat=n=0akk(t),t 0,其中Laguerre函数的展开系数列akk=0满足ak=0eatk(t)dt=(1)k2(a+)k(a)(k+1),k=0,1,.类似地,对于含有一个稳定矩阵A(若A Cnn的所有特征值严格位于复平面的左半平面,则称A为稳定矩阵)的矩阵指数函数eAt,可作下列Laguerre展开16eAt=n=0Akk(t),t 0,第 4 期金辉等：K-power双线性系统基于Laguerre函数的保结构模型降阶方法931其中Laguerre系数矩阵Ak满足Ak=(1)k2(aI+)k(I A)(k+1),k=0,1,I为单位矩阵.3.K-power双线性系统的保结构模型降阶算法考虑如下K-power双线性系统1,17：x1(t)=A1x1(t)+B1u(t),x2(t)=A2x2(t)+mi=1N1ix1(t)ui(t),.xk(t)=Akxk(t)+mi=1Nk1,ixk1(t)ui(t),y(t)=Ckxk(t),(3.1)其中状态变量xj Rnj,Aj Rnjnj,Nji Rnj+1nj,xk Rnk,Ak Rnknk,j=1,2,k 1.输入和输出系数矩阵分别为B1 Rn1m和Ck Rmnk.输出变量y(t)Rm,ui(t)为输入变量u(t)Rm的第i个分量.系统(3.1)的初始条件满足xj(0)=0,j=1,2,k.K-power双线性系统可以转化为具有如下特殊块结构的系统:x1(t)x2(t).xk(t)=A1A2.Akx1(t)x2(t).xk(t)+pi=1000N1i00.0Nk1,i0 x1(t)x2(t).xk(t)ui(t)+B10.0u(t),y(t)=00Ckx1(t)x2(t).xk(t).(3.2)系统(3.2)可以等价转化为如下双线性系统:x(t)=Ax(t)+mi=1Nix(t)ui(t)+Bu(t),y(t)=Cx(t),(3.3)其中x(t)Rn,u(t)Rm,A Rnn,B Rnm,C Rmn,Ni Rnn,n=kj=1nj,i=1,2,m,并且A=A1A2.Ak,Ni=000N1i00.0Nk1,i0,B=B10.0C=00Ck,x(t)=xT1(t)xTk1(t)xTk(t)T.932应用数学2023双线性系统(3.3)的交叉Gram矩阵9为R=i=100PiQTidt1dti,(3.4)其中P1=P1(t1)=eAt1B,Q1=Q1(t1)=eATt1CT,Pi=Pi(t1,.,ti)=eAtiN1Pi1N2Pi1.NmPi1,Qi=Qi(t1,.,ti)=eATtiNT1Qi1NT2Qi1.NTmQi1,其满足如下广义的Sylvester方程:AR+RA+mi=1NiRNi+BC=0.首先,通过Laguerre函数将矩阵eAtB和eATtCT展开为如下近似形式eAtB N1i=0AiBi(t),eATtCTN1i=0ATiCTi(t),其中Ai=(1)i2(I+A)i(I A)(i+1),i=0,1,N 1.因此,可以得到R=0eAtB(eATtCT)Tdt=0(N1i=0AiBi(t)(N1i=0ATiCTi(t)Tdt=0A0BA1BAN1B0(t)I1(t)I.N1(t)I 0(t)I1(t)IN1(t)I(AT0CT)T(AT1CT)T.(ATN1CT)Tdt.根据Laguerre函数的正交性,则有R FGT,其中F=A0BA1BAN1B,G=AT0CTAT1CTATN1CT,且Ai由如下递推公式计算A0=2(I A)1,Ai=(A+I)(A I)1Ai1,i=1,2,N 1.基于以上过程,对双线性系统(3.3)的交叉Gram矩阵R作近似低秩分解,则有R1=0P1QT1dt1=0eAtB(eATtCT)Tdt1 F1GT1,第 4 期金辉等：K-power双线性系统基于Laguerre函数的保结构模型降阶方法933其中F1=F1,0F1,1F1,N1 RnNm,G1=G1,0G1,1G1,N1,且F1,j和G1,j由如下递推公式计算F1,0=2(I A)1B,F1,j=(A+I)(A I)1F1,j1,G1,0=2(I AT)1CT,G1,j=(AT+I)(AT I)1G1,j1,j=1,2,N 1.式(3.4)中级数依赖于双线性系统Volterra级数的核,考虑其前两项级数如下R2=0P1QT1dt1+00P2QT2dt1dt2 F1GT1+00eAt2N1P1N2P1NmP1eATt2NT1Q1NT2Q1.NTmQ1Tdt1dt2=F1GT1+0eAt2N1F1N2F1NmF1eATt2NT1G1NT2G1NTmG1Tdt2 F1GT1+F2GT2,其中F2=F2,0F2,1F2,N1 RnN2m2,G2=G2,0G2,1G2,N1,且F2,j和G2,j由如下递推公式计算F2,0=2(I A)1N1F1N2F1NmF1,F2,j=(A+I)(A I)1F2,j1,G2,0=2(I AT)1NT1G1NT2G1NTmG1,G2,j=(AT+I)(AT I)1G1,j1,j=1,2,N 1.则对于双线性系统的前l项截断的交叉Gram矩阵Rl,得到其如下低秩分解Rl=li=100PiQTidt1dti FGT,(3.5)其中F=F1F2Fl Rnli=1Nimi,Fi=Fi,0Fi,1Fi,N1,G=G1G2Gl,Gi=Gi,0Gi,1Gi,N1,且Fi,j和Gi,j由如下递推公式计算F1,0=2(I A)1B,F1,j=(A+I)(A I)1F1,j1,Fi,0=2(I A)1N1Fi1N2Fi1NmFi1,Fi,j=(A+I)(A I)1Fi,j1,934应用数学2023G1,0=2(I AT)1CT,G1,j=(AT+I)(AT I)1G1,j1,Gi,0=2(I AT)1NT1Gi1NT2Gi1NTmGi1,Gi,j=(AT+I)(AT I)1Gi,j1,i=1,2,l,j=1,2,N 1.进一步,我们给出上述方法与交替方向隐式(ADI)方法18的联系.ADI方法是一种流行的求解Lyapunov和Sylvester方程的低秩Cholesky因子的算法.通过计算可控Gram矩阵P=0eAtBBTeATtdt的低秩因子,P ZjZTj的Cholesky因子ADI(CF-ADI)迭代法如下18z1=21(A+1I)1B,Z1=z1,i=(2i12i)I (i+1+i)(A+i+1I)1,zj=j1zj1,Zj=Zj1zj,j=2,3,N,其中1,2,N为CF-ADI参数,且Re(j)0)为可逆对角矩阵,并且UT1jU1j=VT1jV1j=IrNj,其中IrNj为rNj阶单位矩阵.然后,令Urj和Vrj分别为矩阵U1j和V1j的前rj列.最后,通过构造投影变换矩阵Tj=FjVrj12rj Rnjrj,Sj=12rjUTrjGTj Rrjnj,其中rj是rNj的主子矩阵.第 4 期金辉等：K-power双线性系统基于Laguerre函数的保结构模型降阶方法935可得系统(3.1)的降阶系统如下 x1(t)=A1 x1(t)+B1u(t)x2(t)=A2 x2(t)+mi=1N1i x1(t)ui(t),.xk(t)=Ak xk(t)+mi=1Nk1,i xk1(t)ui(t),y(t)=Ck xk(t),(3.6)其中Aj=SjAjTj,Nji=Sj+1NjiTj,Ak=SkAkTk,B1=S1B1,Ck=CkTk,j=1,2,k1.上述降阶过程可由如下算法1描述：算法1K-power双线性系统基于Laguerre函数的保结构模型降阶算法步0输入:Aj,B1,Ck,Nji,N,l,tolj,n=kj=1nj;步1通过式(3.5)计算低秩分解因子F和G;步2对F和G分块:F=FT1,FT2,FTjT,G=GT1,GT2,GTjT;步3进行SVD分解:GTjFj=UjjVTj,Urj=Uj(:,1:rj),Vrj=Vj(:,1:rj),rj=j(1:rj,1:rj);其中rj由给定的误差限“tolj”自适应选取:=2rNjj=rj+1j tolj;步4计算投影矩阵:Tj=FjVrj12rj,Sj=12rjUTrjGTj;步5构造降阶系统:Aj=SjAjTj,Nji=Sj+1NjiTj,Ak=SkAkTk,B1=S1B1,Ck=CkTk,j=1,2,k 1;步6输出:Aj,Nji,B1,Ck,r=kj=1rj.在实际计算中,N为矩阵指数函数的Laguerre函数展开项数,其一般取值在10到30之间,其最优值还有待于进一步的研究.l是交叉Gram矩阵Rl的截断项,在实际应用中,通常取为2.在算法1中,对于给定的误差限“tolj”,降阶模型的阶数rj可由近似误差估计式=2rNjj=rj+1j tolj自适应地选取,其中1j,2j,rNj为矩阵GTjFj依次由大到小的奇异值.此外,Laguerre参数的取值直接影响降阶系统的精度和稳定性.文13对Laguerre参数的选取作了进一步研究.对于线性系统,参数的最优取值可能在4B,2B区间,其中B为线性系统的带宽,并且文16 给出了Laguerre参数的极大极小优化问题:min0max(A)?+?,(3.7)其中(A)为A的特征值.对于更一般的情况,如何选取最优的Laguerre参数需要进一步研究.数值结果表明,该最优问题对我们算法中参数的选取具有指导作用.降阶系统(3.6)可以写为如下等价的双线性系统形式:x(t)=A x(t)+mi=1Ni x(t)ui(t)+Bu(t),y(t)=C x(t),其中936应用数学2023A=A1A2.Ak,Ni=000N1i00.0Nk1,i0,B=B10.0C=00Ck,x(t)=xT1(t)xTk1(t)xTk(t)T.由K-power双线性系统的等价双线性系统及其降阶系统,可以得到双线性系统和降阶系统的误差系统如下:x(t)=A x(t)+mi=1Ni x(t)ui(t)+Bu(t),y(t)=C x(t),其中 x(t)=xT(t)xT(t)T,A=diagA,A,Ni=diagNi,Ni,B=BTBTT,C=CC并且 y(t)=y(t)y(t).定理3.120对于双线性系统,如果A是稳定的,并且可控Gram矩阵P和可观Gram矩阵Q均存在,则双线性系统的H2范数表示为 H2=tr(CPCT)=tr(BTQB).同理,对于给定的误差系统,如果矩阵A是稳定的,则误差系统的H2范数可以表示为 H2=tr(CPCT)=tr(BTQB),其中,可控Gram矩阵P和可观Gram矩阵Q满足如下Lyapunov方程AP+PAT+mi=1NiPNTi+BBT=0,ATQ+QA+mi=1NTiQNi+CTC=0.对于精确求解双线性系统可控Gram矩阵P和可观Gram矩阵Q的平衡截断方法6和某些H2优化模型降阶方法17,可以得到误差系统的Hankel范数或H2范数.一般而言,对于近似平衡截断方法,即,通过近似求解或分解系统可控Gram矩阵P和可观Gram矩阵Q,很难得到与精确平衡截断方法类似的误差估计式.尽管如此,值得注意的是,XIANG在文21中证明了矩阵指数函数正交多项式近似的指数收敛性,受此启发,我们可以将该结论推广应用到(双)线性系统Gram矩阵的低秩分解的收敛性分析,并由此估计降阶系统与原系统的误差,该问题有待进一步的研究.4.数值实验下面通过一个数值算例来验证上述保结构模型降阶算法的有效性.该数值实验在Intel(R)Core(TM)i5-10210U CPU(1.60 GHZ)PC上进行,内存为16.00GB.同时,选用ode15s函数在MATLAB R2018b中求解系统.考虑由两个子系统组成的单输入单输出(SISO)K-power双线性系统12,17,19,并选取每个第 4 期金辉等：K-power双线性系统基于Laguerre函数的保结构模型降阶方法937子系统的阶数均为3000.该系统的系数矩阵如下:A1=1027102.710,B1=11.1,A2=52252.25,N11=21121.12,C2=0.01T,其中B1,CT2 R3000.通过本文给出的算法1(Algorithm 1)对该系统进行降阶,该系统优化问题(3.7)的解为 10.数值结果表明,该值是有效的.取Laguerre函数展开项数N=4,误差限tol=1010,由此自适应得到的两个降阶子系统的阶数均为4.同时,分别以文22中的Krylov子空间方法(Krylov),文19中的Laguerre正交多项式方法(Laguerre)和文10中的平衡截断方法(BT)对该系统进行降阶.在Krylov子空间方法中,取q1=9匹配该系统第一个子系统的前9阶矩,同时选取p2=2,q2=4匹配第二个子系统的前4阶矩19,由此得到的两个降阶子系统的阶数为9和8.由Laguerre多项式方法和BT方法分别得到的两个降阶子系统的阶数也为9和8.表1分别给出了该系统的输入变量u1(t)=10sin(10t+5)和u2(t)=sin(4t)e0.2t时各模型降阶方法的CPU运行时间(包括降阶的时间和模拟降阶系统的时间)、最大相对误差y(t)yr(t)2/y(t)2和加速比,其中yr(t)为降阶系统的输出响应.表 1各模型降阶方法的CPU运行时间、最大相对误差和加速比方法模型阶数时间(秒)最大相对误差u1(t)u2(t)加速比Original system6000151.021-Algorithm 181.5065.2 1036.4 102100Krylov1717.4141.4 1013.0 1018Laguerre1721.2093.8 1026.4 1017BT17161.8415.4 1035.9 102-图1和图2分别给出了输入变量为u1(t)时各降阶系统和原始系统的瞬态响应及其对应的相对误差.00.511.522.533.544.55t00.10.20.30.40.50.60.70.80.91y(t)Original systemAlgorithm 1KrylovLaguerreBT图1输入为u1时各降阶系统的瞬态响应00.511.522.533.544.55t10-610-510-410-310-210-1100Relative errorsAlgorithm 1KrylovLaguerreBT图2输入为u1(t)时瞬态响应的相对误差图3和图4分别给出了输入变量为u2(t)时各降阶系统和原始系统的瞬态响应及其对应的相对误差.938应用数学202300.511.522.533.544.55t00.010.020.030.040.05y(t)Original systemAlgorithm 1KrylovLaguerreBT图3输入为u2(t)时各降阶系统的瞬态响应00.511.522.533.544.55t10-610-510-410-310-210-1100Relative errorsAlgorithm 1KrylovLaguerreBT图4输入为u2(t)时瞬态响应的相对误差从模拟结果可以看出,对于该SISO K-power双线性系统模型,由算法1得到的降阶系统对原始系统有很好地瞬态响应近似效果.在相同的近似精度下,算法1构造的降阶系统比BT方法构造的降阶系统的阶数小.算法1构造的降阶系统比Krylov子空间方法和Laguerre多项式方法构造的降阶系统有略好的近似精度,并且算法1得到的降阶系统的阶数小于这两种方法(Krylov,Laguerre)得到的降阶系统的阶数.此外,算法1构造和模拟降阶系统的时间少于其他三种方法(Krylov,Laguerre,BT)所用的时间.这些结果表明,本文提出的算法对于该K-power双线性系统是有效的.5.结论本文提出了一种K-power双线性系统基于Laguerre函数的保结构模型降阶方法.该方法首先将K-power双线性系统转化为等价的双线性系统,再结合Laguerre函数的正交性,将矩阵指数函数通过Laguerre函数展开,然后计算出双线性系统交叉Gram矩阵的低秩分解因子,从而构造K-power双线性系统各子系统相应的投影变换,最后得到保结构的降阶模型.通过所提算法得到的降阶系统能很好地近似原始系统的瞬态相应,并且该方法计算高效,具有一定的自适应性.最后,数值实验有效验证了该算法的有效性.参考文献：1 蒋耀林.模型降阶方法M.北京:科学出版社,2010.2 ANTOULAS A C.Approximation of Large-scale Dynamical SystemsM.Philadelphia:SIAM,2005.3 MOORE B.Principal component analysis in linear systems:controllability,observability,and modelreductionJ.IEEE Transactions on Automatic Control,1981,26(1):17-32.4 BAUR U,BENNER P,FENG Lihong.Model order reduction for linear and nonlinear systems:asystem-theoretic perspectiveJ.Archives of Computational Methods in Engineering,2014,21(4):331-358.5 MOHLER R R.Bilinear Control ProcessesM.New York:Academic Press,1973.6 BENNER P,DAMM T.Lyapunov equations,energy functionals,and model order reduction of bilinearand stochastic systemsJ.SIAM journal on control and optimization,2011,49(2):686-711.7 BREITEN T,DAMM T.Krylov subspace methods for model order reduction of bilinear controlsystemsJ.Systems&Control Letters,2010,59(8):443-450.8 XU Kangli,JIANG Yaolin.An unconstrained H2model order reduction optimisation algorithm basedon the Stiefel manifold for bilinear systemsJ.International Journal of Control,2019,92(4):950-959.9 SHAKER H R,TAHAVORI M.Control configuration selection for bilinear systems via generalisedHankel interaction index arrayJ.International Journal of Control,2015,88(1):30-37.10 AL-BAIYAT S A,BETTAYEB M.A new model reduction scheme for K-power bilinear system-sC/Proceedings of 32nd IEEE Conference on Decision and Control.IEEE,1993:22-27.第 4 期金辉等：K-power双线性系统基于Laguerre函数的保结构模型降阶方法93911 XIE Baohui,SYRMOS V L.The-reciprocal system for model reduction of K-power bilinear sys-temsC/Proceedings of 35th IEEE Conference on Decision and Control.IEEE,1996,4:4294-4299.12 WANG Xiaolong,JIANG Yaolin.On model reduction of K-power bilinear systemsJ.InternationalJournal of Systems Science,2014,45(9):1978-1990.13 KNOCKAERT L,DE ZUTTER D.Laguerre-SVD reduced-order modelingJ.IEEE Transactions onMicrowave Theory and Techniques,2000,48(9):1469-1475.14 DAI Liu,XIAO Zhihua,ZHANG Renzheng,et al.Laguerre-Gramian-based structure-preserving mod-el order reduction for second-order form systemsJ.Transactions of the Institute of Measurement andControl,2022,44(5):1113-1124.15 MOORE G.Orthogonal polynomial expansions for the matrix exponentialJ.Linear Algebra and itsApplications,2011,435(3):537-559.16 XIAO Zhihua,SONG Qiuyan,JIANG Yaolin,et al.Model order reduction of linear and bilinearsystems via low-rank Gramian approximationJ.Applied Mathematical Modelling,2022,106:100-113.17 杨平,徐康丽,蒋耀林.K-power系统在Grassmann流形上的单侧模型降阶方法J.数值计算与计算机应用,2016,37(1):41-56.18 LI J R,WHITE J.Low rank solution of Lyapunov equationsJ.SIAM Journal on Matrix Analysisand Applications,2002,24(1):260-280.19 QI Zhenzhong,JIANG Yaolin,XIAO Zhihua.Dimension reduction for K-power bilinear systemsusing orthogonal polynomials and Arnoldi algorithmJ.International Journal of Systems Science,2021,52(10):2048-2063.20 ZHANG Liqian,JAMES L.On H2model reduction of bilinear systemsJ.Automatica,2002,38(2):205-216.21 XIANG Shuhuang.On error bounds for orthogonal polynomial expansions and Gauss-type quadra-tureJ.SIAM Journal on Numerical Analysis,2012,50(3):1240-1263.22 BAI Zhaojun,SKOOGH D.A projection method for model reduction of bilinear dynamical systemsJ.Linear Algebra and its Applications,2006,415(2-3):406-425.Structure Preserving Model Order Reduction for K-powerBilinear Systems Based on Laguerre FunctionsJIN Hui1,XIAO Zhihua1,QI Zhenzhong2(1.School of Information and Mathematics,Yangtze University,Jingzhou434023,China;2.School of Mathematics,Northwest University,Xian 710127,China)Abstract:A structure preserving model order reduction method for K-power bilinear systemsbased on Laguerre functions is proposed.Firstly,the method aims to transform the K-power bilinearsystem into a general bilinear system and calculate the approximate low-rank factors of the cross Gramianof the bilinear system by Laguerre functions expansion of the matrix exponential function.After that,the approximate balanced system of the K-power bilinear system is constructed by the correspondingprojection transformation of each subsystem.Then,the reduced order model is obtained by truncatingthe states corresponding to smaller singular values.The process of model order reduction is flexible,efficient and adaptive.Finally,a numerical experiment is given to demonstrate the effectiveness of thealgorithm.Key words:Model order reduction;K-power bilinear system;Laguerre function;Cross Gramian