将非线性动力系统作为已辨识的线性参数可变模型进行降阶建模
本示例展示了如何使用线性参数可变 (LPV) 建模技术对非线性系统进行降阶建模 (ROM)。用于描述该方法的非线性系统是由多个非线性质量-弹簧-阻尼系统级联而成的。降阶模型的建立基于这样一个思路:通过线性化或线性模型辨识,在非线性模型的工作范围内获得一个线性模型数组。这些局部线性模型通过插值技术被拼接在一起,从而生成一个线性参数可变的表示。此类表示方法既有助于生成仿真速度更快的代理模型,也有助于控制设计。
为了进行局部模型计算,操作空间由一组调度参数(下图中的 p1、p2)来表征。这些参数会根据预先选定的值在网格上进行变化。首先通过模型调整(请参阅 Simulink® Control Design™ 中的 findop)确定与每个网格点对应的平衡运行工况。在每个工作点,都会进行局部扰动试验,以获得系统动态特性的线性逼近。由此获得的局部模型被用于构建基于网格的 LPV 模型。此示例要求所有本地模型具有相同的结构。特别是,它们都必须使用相同的状态定义。本示例采用正交分解 (POD) 方法,将测得的局部响应投影到一个共同的低维空间中。利用预测的状态轨迹,通过动态模式分解 (DMD) 构建低阶状态空间模型。

该示例使用一个非线性质量-弹簧-阻尼器 Simulink 模型作为参考高保真系统。该系统表现出非线性刚度和阻尼力;这些力是质量位移的多项式函数。LPV 模型的开发遵循以下步骤:
该非线性模型是在一系列位移值范围内进行调整的。
在每个调整点,都会进行局部小信号仿真,并记录相关数据。
随后,利用这些数据集构建了一组线性状态空间模型数组。对一个自由度的线性参数可变模型,沿力 的特定轨迹进行了仿真。由此便得到了多自由度降阶模型。
n4sid命令用于根据输入-输出数据估计线性状态空间模型,该命令被用于辨识质量-弹簧-阻尼器系统级联结构中平衡点的线性模型。
质量-弹簧-阻尼器模型
考虑下图所示的单质量-弹簧-阻尼系统 (MSD),其中 为质量, 为阻尼系数, 为弹簧刚度系数, 是质量相对于平衡位置的位移, 和 是作用于质量上的输入力。 代表可测量的扰动, 是可控输入。

MSD 系统可按 [1] 中的模型进行建模:
其中 表示由弹簧刚度产生的力, 表示阻尼力。以位移 和速度 作为状态变量,该方程可表示为状态空间形式:
.
定义 并对上述方程进行离散化,可得
.
通过将 简化为 来简化记法,从而得到:
在此示例中,我们考虑由多项式描述的刚度和阻尼力中的非线性:
,
且阻尼系数由以下公式给出:
,
其中 (加固弹簧)。设 表示用于调度非线性动态特性的时间可变运行工况。例如,在此示例中,您将该行为安排在扰动力的不同值上,即 。如果 和扰动力 保持恒定在 ,则该质量将移动到平衡位置 。可以通过求解以 为变量的三次方程 来确定该平衡位置。请注意,在平衡位置时,阻尼力为零。因此,质量-弹簧-阻尼系统的参数化平衡条件集合为:
.
将该方程以对应于固定值 的平衡条件为基准进行线性化,可得:
其中:
.
在工作条件下 ,线性化的弹簧常数为 。同样,(因为在平衡状态下速度为零)。有关详细信息,请参阅 [1] 中的第 4.2 节。
本示例中用于仿真的模型是由 100 个质量块通过弹簧和阻尼器串联而成的级联系统,如上图所示。
利用扰动试验进行局部分析
设置仿真 MSD 系统所需参数的值。考虑使用:、、、。
% System parameters M = 100; % Number of masses m = 1; % Mass of individual blocks, kg k1 = 0.5; % Linear spring constant, N/m k2 = 1; % Cubic spring constant, N/m^3 b1 = 1; % Linear damping constant, N/(m/s) b2 = 0.1; % Cubic damping constant, N/(m/s)^3 % Model parameters Tf = 20; % simulation time dt = 0.01; % simulation time step (fixed step)
实验设置
如 [1] 所述,在包含 100 个质量点的 MSD 系统上进行了一项扰动试验。首先,为 MSD 系统定义了一组工作点。然后,通过针对每个工作点运行仿真来收集数据。
% Initial condition % (needed to update the size of x in mdl_NL) x0 = zeros(M*2,1); % Simulation times dt = 0.1; % Time step, sec Tf = 20; tin = 0:0.1:Tf; % Time vector for snapshots, sec % Select Parameter Grid F_vec = (0:0.2:2)'; dut = [0 0]; dFt = [0 0]; F = F_vec(1); U0 = 0; Ngrid = numel(F_vec); Nt = numel(tin); % Batch trimming params(1).Name = 'F'; params(1).Value = F_vec;
为 Simulink MSD_NL.slx 模型创建工作点规格。
opspec = operspec('MSD_NL'); opopt = findopOptions('DisplayReport','off'); op = findop('MSD_NL',opspec, params, opopt);
在每个网格点收集数据。
% Collect data at each grid point xtrim = zeros(2*M,Ngrid); dXall = zeros(2*M,Nt,Ngrid); dUall = zeros(1,Nt,Ngrid); dYall = zeros(1,Nt,Ngrid); xOff = zeros(2*M,1,Ngrid); yOff = zeros(1,1,Ngrid); for i=1:Ngrid % Grab current grid point F = F_vec(i); U0 = 0; xtrim(:,i) = op(i).States.x; % Excite nonlinear system with chirp input, u w0 = 0.1; wf = 1+F; win = w0*(wf/w0).^(tin/Tf); uin = 1e-1*sin( win.*tin ); dut = [tin(:) uin(:)]; dFt = [0, 0; Tf, 0]; % Simulate using Simulink model x0 = op(i).States.x; simout_NL_perturb = sim('MSD_NL'); states = simout_NL_perturb.simout.Data; states = (squeeze(states))'; % Collect snapshots dXall(:,:,i) = states'-repmat(x0,[1 Nt]); dUall(:,:,i) = uin; % note that the 100th state is treated as model output dYall(:,:,i) = dXall(M,:,i); % Collect offset data (input offset is 0) xOff(:,:,i) = op(i).States.x; yOff(:,:,i) = op(i).States.x(M); end
降阶线性参数可变建模
使用线性参数可变模型进行降阶建模,其中局部模型是通过动态模态分解 (DMD) 获得的。
X = dXall; Y = dYall; U = dUall; Nx = size(X,1); Nu = size(U,1); Ny = size(Y,1); Ngrid = size(X,3);
将各个网格点的状态轨迹沿时间方向拼接起来,以对状态进行 PCA 型投影。
% Concatenate the Ngrid (=1) state trajectories % along the time (second) dimension. X0all = reshape(X(:,1:end-1,:),Nx,[]); % Plot a few for verification plot(X0all([100 200],:)')

% POD modes using all state snapshots
[UU,SS,~] = svd(X0all);分析 SS 中的奇异值,以确定合适的投影维度。
bar(diag(SS)) xlim([0 10]) ylabel('Singular values') xlabel('Order')

该图大致表明这是一个五阶模型。将网格上获得的状态轨迹投影到 5 维空间中。
% Select modes for subspace projection Q = UU(:,1:5); % Project all the state trajectories into the 5-dimensional space Xs = pagemtimes(Q',X);
预测的状态数据 Xs 以及每个网格点处的相应输入和输出轨迹,可用于辨识算法中以构建局部线性模型。例如,在此示例中,我们将辨识问题表述为对状态和输出轨迹进行一步预测的问题:
这里, 和 是被测量。因此,可以通过线性回归求得未知数 的值。该方法是一种动态模态分解 (DMD) 形式,常用于流体动态特性研究 - DMD 模态和特征值通过振荡组件来描述时间序列中观察到的动态特性。DMD 技术用于在每个网格点上获得独立的线性模型。
% Generate values for Xs(t+1), and Xs(t) over a time span % Time is along the second dimension in the variables Xs, Y, U Xs_future = Xs(:,2:end,:); % Xs(t+1) Xs_current = Xs(:,1:end-1,:); % Xs(t) Y_current = Y(:,1:end-1,:); % Y(t) U_current = U(:,1:end-1,:); % U(t) H = cat(1,Xs_future,Y_current); % response to be predicted R = cat(1,Xs_current,U_current); % regressor data ("predictor") % Compute parameter estimates by linear regression ABCD = pagemrdivide(H,R); A = ABCD(1:5,1:5,:); B = ABCD(1:5,6:end,:); C = ABCD(6:end,1:5,:); D = ABCD(6:end,6:end,:); % Create a state space model array % to represent the dynamics over the entire grid Gred = ss(A,B,C,D,dt)
Gred(:,:,1,1) =
A =
x1 x2 x3 x4 x5
x1 0.9968 -0.1 0.01319 0.009131 -0.004158
x2 0.0598 0.8173 0.2205 0.1001 -0.3144
x3 0.003302 0.0785 0.7367 0.04445 0.2369
x4 -0.00134 0.01687 -0.09067 0.9168 0.09428
x5 -0.0004078 -0.003254 0.04517 0.08203 0.8705
B =
u1
x1 -1.956e-05
x2 0.1129
x3 -0.03799
x4 -0.009784
x5 0.004311
C =
x1 x2 x3 x4 x5
y1 -0.9049 0.08258 0.03128 0.3008 0.2373
D =
u1
y1 -0.0006643
Gred(:,:,2,1) =
A =
x1 x2 x3 x4 x5
x1 0.9969 -0.1001 0.01301 0.008776 -0.004979
x2 0.07983 0.8193 0.2194 0.04261 -0.4338
x3 -0.004858 0.07868 0.7395 0.07153 0.2777
x4 -0.002761 0.01683 -0.08896 0.9178 0.09539
x5 0.0007509 -0.003889 0.04217 0.08096 0.8763
B =
u1
x1 -5.993e-05
x2 0.1087
x3 -0.03629
x4 -0.009364
x5 0.003973
C =
x1 x2 x3 x4 x5
y1 -0.9046 0.08234 0.03654 0.2987 0.2158
D =
u1
y1 -0.0003161
Gred(:,:,3,1) =
A =
x1 x2 x3 x4 x5
x1 0.9972 -0.1003 0.01283 0.008032 -0.00663
x2 0.121 0.8147 0.2297 -0.009268 -0.6053
x3 -0.02094 0.08264 0.7379 0.09681 0.3403
x4 -0.005489 0.01761 -0.08895 0.9193 0.1044
x5 0.002551 -0.005922 0.04063 0.07858 0.8736
B =
u1
x1 -6.534e-05
x2 0.1081
x3 -0.03621
x4 -0.009317
x5 0.004041
C =
x1 x2 x3 x4 x5
y1 -0.9035 0.08259 0.03938 0.2953 0.2003
D =
u1
y1 -0.0001989
Gred(:,:,4,1) =
A =
x1 x2 x3 x4 x5
x1 0.9974 -0.1006 0.01263 0.006901 -0.008381
x2 0.1685 0.8085 0.2354 -0.05807 -0.7645
x3 -0.03888 0.08726 0.7382 0.1219 0.3977
x4 -0.008794 0.01842 -0.08986 0.924 0.1168
x5 0.004541 -0.00794 0.04038 0.07137 0.8663
B =
u1
x1 -7.154e-05
x2 0.1082
x3 -0.03629
x4 -0.009334
x5 0.00413
C =
x1 x2 x3 x4 x5
y1 -0.9023 0.08265 0.04014 0.2957 0.1912
D =
u1
y1 -0.0001015
Gred(:,:,5,1) =
A =
x1 x2 x3 x4 x5
x1 0.9975 -0.1008 0.01234 0.00552 -0.009672
x2 0.2145 0.8042 0.2269 -0.1153 -0.8794
x3 -0.05591 0.09101 0.7432 0.1501 0.438
x4 -0.01218 0.0189 -0.09013 0.9308 0.1265
x5 0.006532 -0.009566 0.0403 0.06192 0.8605
B =
u1
x1 -9.825e-05
x2 0.1067
x3 -0.0358
x4 -0.009218
x5 0.004116
C =
x1 x2 x3 x4 x5
y1 -0.9014 0.08257 0.0395 0.2971 0.1851
D =
u1
y1 -5.823e-05
Gred(:,:,6,1) =
A =
x1 x2 x3 x4 x5
x1 0.9977 -0.1009 0.01201 0.004158 -0.01051
x2 0.2563 0.8014 0.2105 -0.1711 -0.9568
x3 -0.07114 0.09401 0.7507 0.1777 0.4641
x4 -0.01535 0.01915 -0.08987 0.9375 0.1328
x5 0.008356 -0.01088 0.04022 0.05266 0.8574
B =
u1
x1 -0.000138
x2 0.1041
x3 -0.03488
x4 -0.008992
x5 0.004029
C =
x1 x2 x3 x4 x5
y1 -0.9006 0.0824 0.03833 0.2983 0.1799
D =
u1
y1 -3.959e-05
Gred(:,:,7,1) =
A =
x1 x2 x3 x4 x5
x1 0.9978 -0.1011 0.01177 0.002849 -0.01165
x2 0.2972 0.7973 0.1986 -0.2263 -1.06
x3 -0.08582 0.09732 0.7562 0.2044 0.4979
x4 -0.01854 0.01948 -0.09004 0.9438 0.141
x5 0.01012 -0.0121 0.04049 0.04419 0.8545
B =
u1
x1 -0.0001682
x2 0.1021
x3 -0.03416
x4 -0.008821
x5 0.00396
C =
x1 x2 x3 x4 x5
y1 -0.9 0.08218 0.037 0.2987 0.1747
D =
u1
y1 -3.35e-05
Gred(:,:,8,1) =
A =
x1 x2 x3 x4 x5
x1 0.9979 -0.1013 0.01159 0.001579 -0.01285
x2 0.3381 0.7918 0.189 -0.283 -1.173
x3 -0.1003 0.101 0.7606 0.2311 0.5342
x4 -0.02178 0.01991 -0.09041 0.9497 0.1496
x5 0.01183 -0.01327 0.04099 0.03656 0.8523
B =
u1
x1 -0.0001875
x2 0.1008
x3 -0.03371
x4 -0.008723
x5 0.003927
C =
x1 x2 x3 x4 x5
y1 -0.8994 0.08194 0.03562 0.2985 0.1691
D =
u1
y1 -3.29e-05
Gred(:,:,9,1) =
A =
x1 x2 x3 x4 x5
x1 0.9979 -0.1015 0.01134 0.0004775 -0.01326
x2 0.3734 0.7878 0.1718 -0.3315 -1.229
x3 -0.1126 0.1039 0.7674 0.2543 0.5499
x4 -0.02453 0.02018 -0.09004 0.9545 0.1525
x5 0.01323 -0.01427 0.04122 0.03023 0.8538
B =
u1
x1 -0.0002154
x2 0.09881
x3 -0.03301
x4 -0.008556
x5 0.003861
C =
x1 x2 x3 x4 x5
y1 -0.8989 0.08168 0.03425 0.2976 0.1633
D =
u1
y1 -3.545e-05
Gred(:,:,10,1) =
A =
x1 x2 x3 x4 x5
x1 0.9978 -0.1016 0.01102 -0.0004348 -0.01298
x2 0.4018 0.7859 0.1488 -0.3704 -1.237
x3 -0.1225 0.1061 0.776 0.2735 0.5484
x4 -0.02665 0.02024 -0.08908 0.9579 0.1503
x5 0.01423 -0.0151 0.04123 0.02525 0.8588
B =
u1
x1 -0.0002555
x2 0.09582
x3 -0.03197
x4 -0.008294
x5 0.003747
C =
x1 x2 x3 x4 x5
y1 -0.8983 0.08141 0.03287 0.2963 0.1572
D =
u1
y1 -4.082e-05
Gred(:,:,11,1) =
A =
x1 x2 x3 x4 x5
x1 0.9978 -0.1017 0.01076 -0.001259 -0.01285
x2 0.4298 0.783 0.1292 -0.4076 -1.264
x3 -0.132 0.1085 0.7832 0.2914 0.5523
x4 -0.02871 0.02038 -0.08841 0.9608 0.1493
x5 0.01511 -0.0159 0.04148 0.02121 0.8642
B =
u1
x1 -0.0002888
x2 0.09334
x3 -0.03111
x4 -0.008077
x5 0.003654
C =
x1 x2 x3 x4 x5
y1 -0.8978 0.08113 0.03149 0.2945 0.151
D =
u1
y1 -4.865e-05
Sample time: 0.1 seconds
11x1 array of discrete-time state-space models.
Model Properties
% Collect offset information % (obtained earlier during operating point search) Gred.SamplingGrid.F = F_vec; % scheduling parameter grid Gred_offset_x = pagemtimes(Q',xOff); Gred_offset_y = yOff; bodemag(Gred) % frequency responses of the 11 models

LTI 数组 Gred 结合偏移数据可用于建立 LPV 模型。有两种方法可以做到这一点:
在 MATLAB® 中使用
ssInterpolant函数创建一个lvpss模型在 Simulink 中使用 LPV System 模块
LPV 模型的创建与仿真
使用正弦输入力对非线性模型和线性参数可变模型进行仿真,并比较结果。首先,对原始非线性系统进行仿真,以生成参考数据。
Tf = 50; tin = 0:dt:Tf; Nt = numel(tin); % Specify sinusoidal parameter (F) trajectory amp = -1; bias = 1; freq = 0.5; Ft = amp*cos(freq*tin)+bias; dFt = [tin(:) Ft(:)]; F = 0; % Specify input Force dut = [tin(:) zeros(Nt,1)]; dut(tin>=25,2) = 0.1; U0 = 0; % Simulate the high-fidelity nonlinear model simout_NL_MDOF = sim('MSD_NL'); tout = simout_NL_MDOF.tout; simout_NL_MDOF_data = squeeze(simout_NL_MDOF.simout.Data)';
在 Simulink 中对 LPV 模型进行仿真
使用 LPV System 模块来表示 LPV 系统。该模块默认使用线性插值来对线性模型矩阵和偏移量进行插值。
% Simulate ROM LPV open_system('MSD_MDOF_LPV')

simout_NL_MDOF_red = sim('MSD_MDOF_LPV');
simout_NL_MDOF_red_data = squeeze(simout_NL_MDOF_red.simout.Data);
ybar = simout_NL_MDOF_red_data(:,7);比较原始模型 (MSD_NL) 和 LPV 逼近模型 (MSD_MDOF_LPV) 的响应结果。
% Compare the system output (position of the last mass) plot(tout, simout_NL_MDOF_data(:,M),'b',... tout, simout_NL_MDOF_red_data(:,6),'r-.',... tout, ybar,'g:') ylabel('Block 100 Position [m]'); xlabel('Time [sec]'); legend('Nonlinear model',... 'LPV (Simulink)',... 'Output offset (Trim)'); grid on;
![Figure contains an axes object. The axes object with xlabel Time [sec], ylabel Block 100 Position [m] contains 3 objects of type line. These objects represent Nonlinear model, LPV (Simulink), Output offset (Trim).](../../examples/ident/win64/ReducedOrderModelingOfANonlinearMSDSystemUsingIdentExample_07.png)
在 MATLAB 中对 LPV 模型进行仿真
LPV 模型由 MATLAB 中的 lpvss 对象封装。lpvss 对象支持仿真,以及 c2d、feedback 连接等模型操作。当 LPV 模型由一组局部线性模型数组组成时,可以使用 ssInterpolant 命令来创建该 LPV 模型。
xc = squeeze(mat2cell(pagemtimes(Q',xOff),5,1,ones(1,Ngrid))); yc = squeeze(mat2cell(yOff,1,1,ones(1,Ngrid))); Offset = struct('dx',xc,'x',xc,'y',yc,'u',[]); LPVModel = ssInterpolant(Gred, Offset)
Discrete-time state-space LPV model with 1 outputs, 1 inputs, 5 states, and 1 parameters. Model Properties
使用与原始非线性系统相同的输入和调度轨迹来仿真 LPVModel。仿真使用 lsim 命令进行。
Input = dut(:,2)+U0; Time = dut(:,1); Scheduling = dFt(:,2)+F; x0 = zeros(5,1); yLPV = lsim(LPVModel, Input, Time, x0, Scheduling); plot(tout, simout_NL_MDOF_data(:,M),'b',... Time, yLPV,'r-.') ylabel('Block 100 Position [m]'); xlabel('Time [sec]'); legend('Nonlinear model',... 'LPV (MATLAB)'); grid on;
![Figure contains an axes object. The axes object with xlabel Time [sec], ylabel Block 100 Position [m] contains 2 objects of type line. These objects represent Nonlinear model, LPV (MATLAB).](../../examples/ident/win64/ReducedOrderModelingOfANonlinearMSDSystemUsingIdentExample_08.png)
结果表明,原始非线性系统的响应与其通过在调度值网格上进行局部线性建模所获得的 LPV 逼近响应之间具有良好的吻合度。
参考资料
[1] Jennifer Annoni and Peter Seiler."A method to construct reduced‐order parameter‐varying models."International Journal of Robust and Nonlinear Control 27.4 (2017):582-597.