主要内容

本页采用了机器翻译。点击此处可查看英文原文。

将非线性动力系统作为已辨识的线性参数可变模型进行降阶建模

本示例展示了如何使用线性参数可变 (LPV) 建模技术对非线性系统进行降阶建模 (ROM)。用于描述该方法的非线性系统是由多个非线性质量-弹簧-阻尼系统级联而成的。降阶模型的建立基于这样一个思路:通过线性化或线性模型辨识,在非线性模型的工作范围内获得一个线性模型数组。这些局部线性模型通过插值技术被拼接在一起,从而生成一个线性参数可变的表示。此类表示方法既有助于生成仿真速度更快的代理模型,也有助于控制设计。

为了进行局部模型计算,操作空间由一组调度参数(下图中的 p1、p2)来表征。这些参数会根据预先选定的值在网格上进行变化。首先通过模型调整(请参阅 Simulink® Control Design™ 中的 findop)确定与每个网格点对应的平衡运行工况。在每个工作点,都会进行局部扰动试验,以获得系统动态特性的线性逼近。由此获得的局部模型被用于构建基于网格的 LPV 模型。此示例要求所有本地模型具有相同的结构。特别是,它们都必须使用相同的状态定义。本示例采用正交分解 (POD) 方法,将测得的局部响应投影到一个共同的低维空间中。利用预测的状态轨迹,通过动态模式分解 (DMD) 构建低阶状态空间模型。

lpv_rom.png

该示例使用一个非线性质量-弹簧-阻尼器 Simulink 模型作为参考高保真系统。该系统表现出非线性刚度和阻尼力;这些力是质量位移的多项式函数。LPV 模型的开发遵循以下步骤:

  1. 该非线性模型是在一系列位移值范围内进行调整的。

  2. 在每个调整点,都会进行局部小信号仿真,并记录相关数据。

  3. 随后,利用这些数据集构建了一组线性状态空间模型数组。对一个自由度的线性参数可变模型,沿力 F 的特定轨迹进行了仿真。由此便得到了多自由度降阶模型。n4sid 命令用于根据输入-输出数据估计线性状态空间模型,该命令被用于辨识质量-弹簧-阻尼器系统级联结构中平衡点的线性模型。

质量-弹簧-阻尼器模型

考虑下图所示的单质量-弹簧-阻尼系统 (MSD),其中 m 为质量,b 为阻尼系数,k 为弹簧刚度系数,p 是质量相对于平衡位置的位移,Fu 是作用于质量上的输入力。F 代表可测量的扰动,u 是可控输入。

MSD 系统可按 [1] 中的模型进行建模:

mp¨=F+u-k(p)-b(p˙)

其中 k(p) 表示由弹簧刚度产生的力,b(p˙) 表示阻尼力。以位移 p 和速度 v=p˙ 作为状态变量,该方程可表示为状态空间形式:

ddt[pv]=[v1m(F+u-b(v)-k(p))].

定义 x(t)=[p(t)v(t)]T 并对上述方程进行离散化,可得

x(tk+1)=x(tk)+Δt[v(tk)F(tk)+u(tk)-b(v(tk))-k(p(tk))m].

通过将 x(tk) 简化为 xk 来简化记法,从而得到:

xk+1=xk+Δt[vkFk+uk-b(vk)-k(pk)m]

在此示例中,我们考虑由多项式描述的刚度和阻尼力中的非线性:

k(p)=k1p+k2p3,

且阻尼系数由以下公式给出:

b(v)=b1v+b2v3,

其中 k2>0(加固弹簧)。设 ρ(tk)=ρk 表示用于调度非线性动态特性的时间可变运行工况。例如,在此示例中,您将该行为安排在扰动力的不同值上,即 ρk=Fk。如果 u0 和扰动力 F 保持恒定在 Fρ,则该质量将移动到平衡位置 p。可以通过求解以 p 为变量的三次方程 ρ=k(p) 来确定该平衡位置。请注意,在平衡位置时,阻尼力为零。因此,质量-弹簧-阻尼系统的参数化平衡条件集合为:

(x(ρ),u(ρ))([p0],0).

将该方程以对应于固定值 ρ 的平衡条件为基准进行线性化,可得:

xk+1xk+1+A(ρk(xk-xk)+B(ρk)(uk-uk)

其中:

xk=x(ρk),uk=u(ρk),A(ρ)=I+Δt[01-klin(ρ)/m-blin/m],B(ρ)=Δt[01/m].

在工作条件下 ρ,线性化的弹簧常数为 klin(ρ)=k1+3k2p(ρ)2。同样,blin(ρ)=b1(因为在平衡状态下速度为零)。有关详细信息,请参阅 [1] 中的第 4.2 节。

本示例中用于仿真的模型是由 100 个质量块通过弹簧和阻尼器串联而成的级联系统,如上图所示。

利用扰动试验进行局部分析

设置仿真 MSD 系统所需参数的值。考虑使用:m=1kgk1=0.5N/mk2=1N/m3b1=1N/(m/s),b2=0.1N/(m/s)3

% 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],:)')

Figure contains an axes object. The axes object contains 2 objects of type line.

% POD modes using all state snapshots
[UU,SS,~] = svd(X0all);

分析 SS 中的奇异值,以确定合适的投影维度。

bar(diag(SS))
xlim([0 10])
ylabel('Singular values')
xlabel('Order')

Figure contains an axes object. The axes object with xlabel Order, ylabel Singular values contains an object of type bar.

该图大致表明这是一个五阶模型。将网格上获得的状态轨迹投影到 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 以及每个网格点处的相应输入和输出轨迹,可用于辨识算法中以构建局部线性模型。例如,在此示例中,我们将辨识问题表述为对状态和输出轨迹进行一步预测的问题:

Xs(t+1)=AXs(t)+BU(t)y(t)=CXs(t)+DU(t)

这里,Xs(t)U(t) 是被测量。因此,可以通过线性回归求得未知数 A,B,C,D 的值。该方法是一种动态模态分解 (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

MATLAB figure

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')

lpvblk.png

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).

在 MATLAB 中对 LPV 模型进行仿真

LPV 模型由 MATLAB 中的 lpvss 对象封装。lpvss 对象支持仿真,以及 c2dfeedback 连接等模型操作。当 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).

结果表明,原始非线性系统的响应与其通过在调度值网格上进行局部线性建模所获得的 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.