CN104614985B - 一种基于非线性规划的高阶系统最优降阶方法 - Google Patents

Abstract

一种基于非线性规划的高阶系统最优降阶方法，该方法有六大步骤：步骤1：降阶问题描述；步骤2：系统最优降阶指标；步骤3：Arnoldi降阶方法；步骤4：非线性规划求解最优插值点；步骤5：仿真实验检验降阶性能；步骤6：设计结束。第一步确定了高阶线性系统降阶的目的和数学描述；第二步确立了系统降阶的误差范数指标，为最优降阶方法的提出做准备；第三步得到了基于krylov子空间的Arnoldi降阶方法；第四步提出了基于非线性规划的最优插值点求解方法，得到最优降阶模型；第五步是对设计的系统最优降阶方法进行仿真验证；经过上述各步骤后，设计结束。本发明用于处理高阶线性系统，简化控制律设计。

Description

一种基于非线性规划的高阶系统最优降阶方法

技术领域

本发明涉及一种基于非线性规划的高阶系统最优降阶方法，它是一种线性时不变系统的最优降阶方法，是针对单输入单输出系统给出的一种使得降阶前后系统频域响应误差最小的模型降阶方法，用于处理高阶线性系统，简化控制律设计，属于自动控制技术领域。

背景技术

伴随着现代工程技术的飞速发展，出现了诸如飞行控制系统、电力系统、超大规模集成电路等复杂的高阶系统。描述这些线性系统的微分方程个数众多，给系统的数值仿真和控制设计带来巨大的挑战。通过对这些高阶系统进行有效的模型降阶处理可以降低系统分析的难度，减少计算负载，方便仿真模拟。降阶是指通过较少个数的微分方程描述的低阶系统，在特定的频率区间上近似高阶系统，保证降阶前后两个系统的动态响应尽可能接近。

系统降阶的关键在于如何找到一个简单的低阶系统来逼近复杂的高阶系统。两个系统的逼近程度可以通过二者在频域上的传递函数脉冲响应偏差来衡量。现有的模型降方法包括Pade近似法、Routh近似法、奇异值分解的平衡法、Krylov子空间模型降阶方法。其中前两种方法基于系统的暂态响应或者稳态响应，大大限制了其应用范围。奇异值分解的平衡法根据原系统的奇异值进行截断降解，能够确保降阶后系统的稳定性，但是在实施过程中需要求解两个高阶次的Lyapunov方程。Krylov子空间方法是最常用的降阶方法，它通过构造适当的变换矩阵，使得变换矩阵的列向量所张成的空间能包含在适当的Krylov子空间中，保证得到的降阶系统传递函数能与原系统传递函数在特定插值点处的项匹配。由于存在无穷多个可行的插值点，在不同的插值点处采用同样的Krylov子空间降阶方法得到的模型是不一样的，相应的系统响应品质和原来的系统差别也是很大的。可以通过非线性规划方法，求解最优插值点，找到系统的最优降阶模型。

发明内容

1、发明目的

本发明的目的是提供一种基于非线性规划的高阶系统最优降阶方法，它是：针对单输入单输出的线性时不变系统，构建系统降解的误差指标，利用非线性规划策略求解满足最优指标的降阶插值点，最后得到最优降阶模型的传递函数。

2、技术方案

为了达到上述目的，本发明结合流程框图1中的步骤，具体介绍该设计方法的技术方案。

本发明是一种基于非线性规划的高阶系统最优降阶方法，该方法具体步骤如下：

步骤1：降阶问题描述

本发明所针对的是高阶单输入单输出线性系统。高维系统的分析和控制律设计都十分困难。通过适当的降阶方法，使得低阶模型在相同的脉冲输入激励情况下，输出响应尽可能接近原来系统，在一定频域范围内代替原来系统。

步骤2：系统最优降阶指标

尽管系统降阶的方法各不相同，但是可以定义一个通用的范数误差指标来衡量降阶前后系统传递函数输出品质差异。对比分析范数误差的各个组成部分，要求两个系统在降阶模型传递函数G_m(s)的极点镜像对称处即点处越接近越好。

步骤3：Arnoldi降阶方法

基于krylov子空间的Arnoldi降阶方法，可以在不计算传递函数展开式各项系数的情况下，满足项匹配的要求。给出原系统的状态矩阵A,B,C，插值点σ＝{σ₁,σ₂,…,σ_m}，降阶误差阈值ε＝10-⁸，初始矩阵V_j＝[]，初始下标变量j＝0。经过降阶算法，先构造krylov投影子空间V，之后得到降阶模型A_m＝V^TAV,B_m＝V^TB,C_m＝CV。式中符号说明如下：A_m代表降阶以后的状态矩阵，B_m代表降阶以后的输入矩阵，C_m代表降阶以后的输出矩阵。

步骤4：非线性规划求解最优插值点

对于不同的插值点，采用Arnoldi降阶算法得到的降阶模型是不一样的，因此对于原来系统的逼近程度也是不一样的。如果采用传统的牛顿迭代方法，只能求解得到一个插值点。但是通过非线性规划函数Fmincon可以找到一个插值点集合，满足系统最优降阶指标所定义的三个必要条件。求出插值点对应的特征方程系数向量{P₁,P₂,…,P_m}，然后利用Matlab自带的Roots()函数求解对应的根即为插值点{σ₁,σ₂,…,σ_m}。然后采用Arnoldi降阶方法，求出所有插值点对应的降阶模型。使得二范数指标η＝||G(s)-G_m(s)||₂最小的插值点，为需要求解的最优插值点。最优插值点对应的降阶模型，即为系统最优降阶模型。式中符号说明如下：G(s)和G_m(s)分别代表原来高阶系统的传递函数和降阶以后系统的传递函数，||·||₂表示·的二范数。

步骤5：仿真实验检验降阶性能

为了验证所提出的系统最优降阶方法的有效性，检验降阶后系统的响应品质和原模型的偏差，利用matlab7.0仿真软件，进行仿真实验。首先选取一个4阶随机模型，将其降低为1阶，分析发现找到了多个插值点，克服了牛顿迭代方法只能找到一个插值点的缺点，并且在低频和高频段都有较好的逼近效果。之后选取一个标准的测试模型，Build PDE。原来高阶系统n＝120，降阶为m＝6的低阶系统。

步骤6：设计结束

整个设计过程分为六大步骤。第一步确定了高阶线性系统降阶的目的和数学描述；第二步确立了系统降阶的误差范数指标，为最优降阶方法的提出做准备；第三步得到了基于krylov子空间的Arnoldi降阶方法；第四步提出了基于非线性规划的最优插值点求解方法，得到最优降阶模型；第五步是对设计的系统最优降阶方法进行仿真验证。经过上述各步骤后，设计结束。

3、优点及功效

本发明的优点在于针对工程领域常见的高阶数学模型，提供了一种系统最优降阶方法，可以找到一个在频域上尽可能逼近原系统响应品质的低阶系统，显著地降低系统数值仿真和控制律设计难度。通过非线性规划方法，找到了最优降阶模型，克服了传统的牛顿迭代方法只能找到单一插值点的缺点，在高频和低频段都能达到很好的降阶效果。

附图说明

图1本发明实施步骤流程框图。

图2本发明中非线性规划求解插值点原理示意图。

图3本发明随机系统降阶仿真示意图。

图4本发明标准Build PDE系统降阶仿真示意图。

图中符号说如下：

图2中P₁,P₂,P₃分别表示约束方程式(17)的三个不同根。

图3中G(s),G₁(s),G₂(s),G₃(s)分别表示原来4阶随机系统以及3个降阶系统。

图4中G₀(s),G_0i(s)(i＝1,…,5)分别原来标准Build PDE系统系统以及5个降阶系统。

具体实施方式

下面将结合技术方案和附图对本发明做进一步的详细说明。见图1，本发明是一种基于非线性规划的高阶系统最优降阶方法，该方法具体步骤如下：

步骤1：系统降阶问题描述

考虑如下单输入单输出线性时不变系统：

其中为控制输入。系统(1)对应的传递函数为

G(s)＝C(sI_n-A)^-1B (2)

其中，I_n为n阶单位矩阵。假设降阶过后的系统模型为Σ_m，

与之对应的传递函数为：

G_m(s)＝C_m(sI_m-A_m)^-1B_m (4)

这里I_m为m阶单位矩阵。对于系统Σ_m，m≤n,控制输入降阶的目的就是使得Σ和Σ_m在相同的脉冲输入激励情况下，输出响应y(t)和y_m(t)尽可能接近。

步骤2：系统最优降阶指标

首先假设传递函数G(s)和G_m(s)的极点(即A和A_m矩阵的特征值)分别为λ_i和令φ_i,分别表示G(s)和G_m(s)在该点的留数，即有φ_i＝res[G(s),λ_i](i＝1,…,n)，则有

式中符号说明如下：表示二范数的平方，和分别表示从第一项到第m项、从第一项到第n项求和。

以上的二范数指标是一个通用的降阶频域数值指标，不论降阶的方法如何。由于上式中第一项的值很小，G(s)和G_m(s)H₂指标值误差指标主要取决于上式第2项，即因此，降阶前后的系统需要在处越接近越好。已知系统G(s)，σ＝{σ₁,…,σ_m}表示插值点的集合，用这些插值点获得降阶模型G_m(s)，G_m(s)在插值点σ＝{σ₁,…,σ_m}满足项匹配G(s)和一阶导数G'(s)。存在一组插值点σ_opt＝{σ_opt1,…,σ_optm}，满足以下3个必要条件，使得最小。

(a)G_m(σ_opti)＝G(σ_opti)(i＝1,…,m)；

(b)G'_m(σ_opti)＝G'(σ_opti)(i＝1,…,m)；

(c)λ(σ_opt)＝-σ_opt。

这里λ(σ_opt)表示最优降阶模型的极点，也就是最优降阶模型的传递函数特征值。

步骤3：Arnoldi降阶方法

将原来的传递函数进行泰勒级数展开，可得同理可以得到若要达到G(s)与G_m(s)的展开前l+1项匹配的要求，此时必须使得两个系统的马尔科夫参数相等，即CA^jB＝C_mA_m ^jB_m(j＝0,1,2,…,l)。

将G(s)在复平面的某点进行泰勒展开，可得

进一步可以得到

其中为系统Σ在点σ展开式的第j项系数，

由G(s)＝C(sI_n-A)^-1B，可得因此

同理可得

为了达到在插值点传递函数展开式项匹配的目的，降阶以后的系统Σ_m必须满足

的前l+1项和原系统式(6)一样，即泰勒级数展开项的系数满足

式(9)即为系统降阶的项匹配条件。将上式写成响应的状态空间表达形式

同理，对于Σ_m有

因此系统降阶的项匹配条件写成响应的状态空间表达式为

C(σI_n-A)^-(j+1)B＝C_m(σI_m-A_m)^-(j+1)B_m (13)

其中j＝0,1,2,…,l。

显然如果通过逐项计算来达到项匹配的条件，将会加重计算负载，而且这是个后验问题。在不知道降阶模型的情况下，(12)式的右边不能计算出来。基于Krylov子空间的有理Arnoldi方法，可以在不计算各项系数具体值的情况下，通过构造特殊的投影子空间，达到项匹配的目的。首先假设降阶插值点σ＝{σ₁,σ₂,…,σ_m}构造投影子空间V

其中κ_i(A,B；σ_i)＝[B,(σ_iI_n-A)^-1B]，i＝1,2,…,m。

那么降阶以后的系统为

A_m＝V^TAV,B_m＝V^TB,C_m＝CV (15)

则系统Σ_m满足A_m＝V^TAV,B_m＝V^TB,C_m＝CV，那么对于i＝0,1，满足降阶条件C(sI_n-A)^-(i+1)B＝C_m(sI_m-A_m)^-(i+1)B_m，即系统Σ_m在σ处匹配前两项。有理Arnoldi降阶算法流程如下：

3-1.初始化原模型的信息A,B,C，插值点σ＝{σ₁,σ₂,…,σ_m}

降阶误差阈值ε＝10^-5，参量V_j＝[],j＝0。

3-2.构造krylov子空间V

for i＝1→m

ifσ(i)＝∞

else

end

ifStop

end

V_j+1＝[V_j v_j+1]

j＝j+1

ifσ(i)＝∞

else

end

ifStop

end

V_j+1＝[V_j v_j+1]

j＝j+1

end

3-3.输出V

3-4.根据V得到降阶模型Σ_m。A_m＝V^TAV,B_m＝V^TB,C_m＝CV。

步骤4：非线性规划求解最优插值点

通过步骤3中最优降阶插值点的第三个必要条件，发现降阶前后模型在降阶模型的特征值相反数处越接近越好。根据柯西展开定理，最优插值点S是如下方程的根s^m-p₁s^m-1+p₂s^m-2-…+(-1)^m-1p_m-1s+(-1)^mp_m＝0，即

为得到最优插值点S，首先需要求出方程组系数向量P＝(p₀,p₁,…,p_m)，p₀＝1。系数向量的约束方程f(P)为

其中

非线性方程组(16)是由m个变量，m个非线性方组构成的，存在多组解{P₁,P₂,…,P_m}。每一个单独的系数解都对应一个插值点S_i。单纯使用牛顿迭代方法求f(P)＝0具有局限性，易陷入局部最优。每次牛顿迭代随机选取的初值有可能都在某一区域内，导致虽经多次运行，仍无法求出其它区域内的解。如图2所示，有可能每次牛顿迭代随机选取的初值都在P₁所在的那个圆内，则可能无法求得另两组解P₂和P₃。采用Fmincon函数进行非线性规划求解，首先根据非线性约束方程求出范数最小的P₁，然后利用Fmincon求出剩余的m-1个解，也就是{P₂,…,P_m}。有约束非线性最优化问题的一般描述为其中x＝[x₁,x₂,…,x_n]^T，该数学表达式的涵义是在约束G(x)≤0的情况下，求取一组x向量，使函数F(x)取得最小值。值得一提的是，这里的约束既可是等式约束也可是不等式约束。约束条件还可以进一步简化为线性不等式约束A_eqx＝B_eq，线性不等式约束Ax≤B，x变量的上界向量x_M和下界向量x_m，使得x_m≤x≤x_M。Matlab优化工具箱为该函数的调用提供了标准格式为，

[x,fopt,key,c]＝fmincon(F,x0,A,B,Aeq,Beq,x_m,x_M,CFun,OPT) (18)

其中F为目标函数的M文件函数，x0为初始搜索点，各个矩阵约束如果不存在，则用空位矩阵来占位。CFun为给非线性约束函数的M文件函数，OPT为控制选项。最优运算完成后，结果将在变量x中返回，最优化的目标函数将在fopt变量中返回。采用Fmincon函数，主要通过目标函数和约束函数的设计，来找到满足条件的P值。基于Fmincon函数的非线性规划求解最优插值点算法步骤为

4-1.产生随机值P₀，带有约束非线性函数为采用Fmincon非线性规划方法求出P₁；

4-2.带有约束非线性函数为采用Fmincon非线性规划函数求P_i；

4-3.根据式16，已知方程系数向量P_i，采用ROOTS()函数求解对应的插值点σ_i。

4-4.采用最优插值点，根据步骤3中的Arnoldi降阶方法求出所有的降阶模型，然后利用指标min(η＝||G-G_m(s)||₂)，得到最优降阶模型。

步骤5：仿真实验检验降阶性能

仿真时，首先选择随机模型进行降阶，测试所设计算法的降阶性能。原系统的阶数为n＝4，模型信息矩阵A，B，C分别为，

C＝[0.7533 0.0123 0.0491 -0.5085]

对应的传递函数为降阶以后系统的阶数为1，采用步骤4的非线性规划求解最优插值点方法，计算出来的插值点系数分别为{P₁,P₂,P₃}＝{0.9427,10000,-1.54*10^-4}，采用ROOTS求解出来的插值点为{σ₁,σ₂,σ₃}＝{0.9427,10000,-1.54*10^-4}。采用步骤3中的Arnoldi降阶方法，降阶传递函数分别为：和原函数的二范数误差分别为η₁＝0.122，3η₂＝0.1619，η₃＝0.1227。最后得到的最优降阶模型为第一个：原随机系统和3个降阶系统的频域响应伯德图如图3所示。

再选取标准的系统降阶测试模型Build PDE model，原始模型阶数为n＝120，传递函数为G₀(s)降阶为m＝6，降阶后的模型传递函数分别为G₀₁(s)到G₀₅(s)，分析结果如表5-1。本发明标准Build PDE系统和5个降阶系统的频域响应伯德图如图4所示。可以看出，在较大的系统阶数降低情况下，所设计的算法也能够达到对原来高阶系统的频域响应逼近，得到的降阶系统可以在不同频域范围内代替原来高阶系统进行数值模拟和控制算法设计。

表5-1 Build PDE model模型降阶结果分析

骤6：设计结束

整个设计过程，从工程领域的大规模高阶系统入手。首先定义了系统降阶问题的数学描述，给出了评价系统降阶品质的评价指标。利用基于krylov子空间的Arnoldi降阶方法，可以根据插值点求出对应的降阶模型。提出了基于非线性规划Fmincon方法，求解最优插值点，得到系统最优降阶模型。

Claims (1)

1.一种基于非线性规划的高阶系统最优降阶方法，其特征在于：该方法具体步骤如下：

步骤1：降阶问题描述

这里所针对的是高阶单输入单输出线性系统，高维系统的分析和控制律设计都十分困难，通过降阶的方法，使得低阶模型在相同的脉冲输入激励情况下，输出响应尽可能接近原来系统，并在预定频域范围内代替原来系统；

步骤2：系统最优降阶指标

尽管系统降阶的方法各不相同，但是定义一个通用的范数误差指标来衡量降阶前后系统传递函数输出品质差异；对比分析范数误差的各个组成部分，要求两个系统在降阶模型传递函数G_m(s)的极点镜像对称处即点处越接近越好；

步骤3：Arnoldi降阶方法

基于krylov子空间的Arnoldi降阶方法，在不计算传递函数展开式各项系数的情况下，满足项匹配的要求；给出原系统的状态矩阵A，输入矩阵B，输出矩阵C，插值点σ＝{σ₁,σ₂,…,σ_m}，降阶误差阈值ε＝10^-8，初始矩阵V_j＝[]，初始下标变量j＝0；经过降阶算法，先构造krylov投影子空间V，之后得到降阶模型A_m＝V^TAV,B_m＝V^TB,C_m＝CV；式中符号说明如下：A_m代表降阶以后的状态矩阵，B_m代表降阶以后的输入矩阵，C_m代表降阶以后的输出矩阵；

步骤4：非线性规划求解最优插值点

对于不同的插值点，采用Arnoldi降阶算法得到的降阶模型是不一样的，因此对于原来系统的逼近程度也是不一样的，如果采用传统的牛顿迭代方法，只能求解得到一个插值点；但是通过非线性规划函数Fmincon能找到一个插值点集合，满足系统最优降阶指标所定义的三个必要条件；求出插值点对应的特征方程系数向量{P₁,P₂,…,P_m}，然后利用Matlab自带的Roots()函数求解对应的根即为插值点{σ₁,σ₂,…,σ_m}，然后采用Arnoldi降阶方法，求出所有插值点对应的降阶模型，使得二范数指标η＝||G(s)-G_m(s)||₂最小的插值点，为需要求解的最优插值点；最优插值点对应的降阶模型，即为系统最优降阶模型；式中符号说明如下：G(s)和G_m(s)分别代表原来高阶系统的传递函数和降阶以后系统的传递函数，||·||₂表示·的二范数；

步骤5：仿真实验检验降阶性能

为了验证所提出的系统最优降阶方法的有效性，检验降阶后系统的响应品质和原模型的偏差，利用matlab7.0仿真软件，进行仿真实验；首先选取一个4阶随机模型，将其降低为1阶，分析发现找到了多个插值点，克服了牛顿迭代方法只能找到一个插值点的缺点，并且在低频和高频段都有较好的逼近效果，之后选取一个标准的测试模型，Build PDE，原来高阶系统n＝120，降阶为m＝6的低阶系统；

步骤6：设计结束

整个设计过程分为六大步骤：第一步确定了高阶线性系统降阶的目的和数学描述；第二步确立了系统降阶的误差范数指标，为最优降阶方法的提出做准备；第三步得到了基于krylov子空间的Arnoldi降阶方法；第四步提出了基于非线性规划的最优插值点求解方法，得到最优降阶模型；第五步是对设计的系统最优降阶方法进行仿真验证；经过上述各步骤后，设计结束。