空气动力学飞行体建模
本示例展示了对一个大型且复杂的非线性系统进行灰盒建模的过程。其目的是展示利用 IDNLGREY 模型对一个具有众多输入 (10) 和输出 (5) 的系统中的大量参数 (16) 进行估计的能力。该系统是一个空气动力学飞行体。我们建立了一个模型，该模型利用物体的速度（平动速度和角速度）测量值以及与其控制面相关的各种角度，来预测该物体的加速度和速度。
输入-输出数据
我们从名为 aerodata.mat 的数据文件中读取测得的速率、角度和动压：
load('aerodata');该文件包含一个均匀采样的数据集，其中变量 u 和 y 各有 501 个样本。采样时间为 0.02 秒。该数据集是根据另一个更复杂的空气动力学飞行体模型生成的。
接下来，我们创建一个 IDDATA 对象来表示和存储数据：
z = iddata(y, u, 0.02, 'Name', 'Data'); z.InputName = {'Aileron angle' 'Elevator angle' ... 'Rudder angle' 'Dynamic pressure' ... 'Velocity' ... 'Measured angular velocity around x-axis' ... 'Measured angular velocity around y-axis' ... 'Measured angular velocity around z-axis' ... 'Measured angle of attack' ... 'Measured angle of sideslip'}; z.InputUnit = {'rad' 'rad' 'rad' 'kg/(m*s^2)' 'm/s' ... 'rad/s' 'rad/s' 'rad/s' 'rad' 'rad'}; z.OutputName = {'V(x)' ... % Angular velocity around x-axis 'V(y)' ... % Angular velocity around y-axis 'V(z)' ... % Angular velocity around z-axis 'Accel.(y)' ... % Acceleration in y-direction 'Accel.(z)' ... % Acceleration in z-direction }; z.OutputUnit = {'rad/s' 'rad/s' 'rad/s' 'm/s^2' 'm/s^2'}; z.Tstart = 0; z.TimeUnit = 's';
查看输入数据：
figure('Name', [z.Name ': input data'],... 'DefaultAxesTitleFontSizeMultiplier',1,... 'DefaultAxesTitleFontWeight','normal',... 'Position',[50, 50, 850, 620]); for i = 1:size(z.InputData,2) subplot(size(z.InputData,2)/2, 2, i); plot(z.SamplingInstants, z.InputData(:,i)); title(['Input #' num2str(i) ': ' z.InputName{i}],'FontWeight','normal'); xlabel(''); axis tight; if (i > size(z.InputData,2)-2) xlabel([z.Domain ' (' z.TimeUnit ')']); end end
图 1：输入信号。
查看输出数据：
figure('Name', [z.Name ': output data']); h_gcf = gcf; Pos = h_gcf.Position; h_gcf.Position = [Pos(1), Pos(2)-Pos(4)/2, Pos(3), Pos(4)*1.5]; for i = 1:size(z.OutputData,2) subplot(size(z.OutputData,2), 1, i); plot(z.SamplingInstants, z.OutputData(:,i)); title(['Output #' num2str(i) ': ' z.OutputName{i}]); xlabel(''); axis tight; end xlabel([z.Domain ' (' z.TimeUnit ')']);
图 2：输出信号。
乍一看，在输入向量中包含某些输出结果的测量值，这似乎有些奇怪。然而，用于生成数据的模型包含多个积分器，这往往会导致仿真行为不稳定。为避免这种情况，将部分输出信号的测量值通过非线性观测器进行反馈。这些是数据集 z 中的第 6 至第 8 个输入。因此，这是一个闭环系统，建模工作的目标是利用当前和过去行为的测量数据来预测这些输出未来的取值。
系统建模
为了对感兴趣的动态特性建模，我们使用一个 IDNLGREY 模型对象来表示该系统的状态空间结构。通过运用牛顿的基本力学定律和动量定律（平衡方程），可以得到一个合理的结构。为了全面描述模型结构，还采用了基本的气体动力学关系（本构关系）。
C MEX 文件 aero_c.c 通过状态方程、输出方程和初始条件来描述该系统，具体内容如下所述。在此省略运动方程的推导细节，仅展示最终状态方程和输出方程，并注意到其结构相当复杂且具有非线性。
/* State equations. */
void compute_dx(double *dx, double *x, double *u, double **p)
{
/* Retrieve model parameters. */
double *F, *M, *C, *d, *A, *I, *m, *K;
F = p[0]; /* Aerodynamic force coefficient. */
M = p[1]; /* Aerodynamic momentum coefficient. */
C = p[2]; /* Aerodynamic compensation factor. */
d = p[3]; /* Body diameter. */
A = p[4]; /* Body reference area. */
I = p[5]; /* Moment of inertia, x-y-z. */
m = p[6]; /* Mass. */
K = p[7]; /* Feedback gain. */
/* x[0]: Angular velocity around x-axis. */
/* x[1]: Angular velocity around y-axis. */
/* x[2]: Angular velocity around z-axis. */
/* x[3]: Angle of attack. */
/* x[4]: Angle of sideslip. */
dx[0] = 1/I[0]*(d[0]*A[0]*(M[0]*x[4]+0.5*M[1]*d[0]*x[0]/u[4]+M[2]*u[0])*u[3]-(I[2]-I[1])*x[1]*x[2])+K[0]*(u[5]-x[0]);
dx[1] = 1/I[1]*(d[0]*A[0]*(M[3]*x[3]+0.5*M[4]*d[0]*x[1]/u[4]+M[5]*u[1])*u[3]-(I[0]-I[2])*x[0]*x[2])+K[0]*(u[6]-x[1]);
dx[2] = 1/I[2]*(d[0]*A[0]*(M[6]*x[4]+M[7]*x[3]*x[4]+0.5*M[8]*d[0]*x[2]/u[4]+M[9]*u[0]+M[10]*u[2])*u[3]-(I[1]-I[0])*x[0]*x[1])+K[0]*(u[7]-x[2]);
dx[3] = (-A[0]*u[3]*(F[2]*x[3]+F[3]*u[1]))/(m[0]*u[4])-x[0]*x[4]+x[1]+K[0]*(u[8]/u[4]-x[3])+C[0]*pow(x[4],2);
dx[4] = (-A[0]*u[3]*(F[0]*x[4]+F[1]*u[2]))/(m[0]*u[4])-x[2]+x[0]*x[3]+K[0]*(u[9]/u[4]-x[4]);
}
/* Output equations. */
void compute_y(double *y, double *x, double *u, double **p)
{
/* Retrieve model parameters. */
double *F, *A, *m;
F = p[0]; /* Aerodynamic force coefficient. */
A = p[4]; /* Body reference area. */
m = p[6]; /* Mass. */
/* y[0]: Angular velocity around x-axis. */
/* y[1]: Angular velocity around y-axis. */
/* y[2]: Angular velocity around z-axis. */
/* y[3]: Acceleration in y-direction. */
/* y[4]: Acceleration in z-direction. */
y[0] = x[0];
y[1] = x[1];
y[2] = x[2];
y[3] = -A[0]*u[3]*(F[0]*x[4]+F[1]*u[2])/m[0];
y[4] = -A[0]*u[3]*(F[2]*x[3]+F[3]*u[1])/m[0];
}
我们还必须为这 23 个参数提供初始值。我们将部分参数（空气动力学力系数、空气动力学动量系数和惯性矩因子）作为向量，分别定义在 8 个不同的参数对象中。初始参数值部分是通过物理推导得出的，部分则是通过定量推测得出的。最后 4 个参数（A、I、m 和 K）基本上是常数，因此通过固定这些参数，我们得到一个包含 16 个自由参数的模型结构，这些参数分布在参数对象 F、M、C 和 d 中。
Parameters = struct('Name', ... {'Aerodynamic force coefficient' ... % F, 4-by-1 vector. 'Aerodynamic momentum coefficient' ... % M, 11-by-1 vector. 'Aerodynamic compensation factor' ... % C, scalar. 'Body diameter' ... % d, scalar. 'Body reference area' ... % A, scalar. 'Moment of inertia, x-y-z' ... % I, 3-by-1 vector. 'Mass' ... % m, scalar. 'Feedback gain'}, ... % K, scalar. 'Unit', ... {'1/(rad*m^2), 1/(rad*m^2), 1/(rad*m^2), 1/(rad*m^2)' ... ['1/rad, 1/(s*rad), 1/rad, 1/rad, ' ... '1/(s*rad), 1/rad, 1/rad, 1/rad^2, ' ... '1/(s*rad), 1/rad, 1/rad'] ... '1/(s*rad)' 'm' 'm^2' ... 'kg*m^2, kg*m^2,kg*m^2' 'kg' '-'}, ... 'Value', ... {[20.0; -6.0; 35.0; 13.0] ... [-1.0; 15; 3.0; -16.0; -1800; -50; 23.0; -200; -2000; -17.0; -50.0] ... -5.0, 0.17, 0.0227 ... [0.5; 110; 110] 107 6}, ... 'Minimum',... {-Inf(4, 1) -Inf(11, 1) -Inf -Inf -Inf -Inf(3, 1) -Inf -Inf}, ... % Ignore constraints. 'Maximum', ... {Inf(4, 1) Inf(11, 1) Inf Inf Inf Inf(3, 1) Inf Inf}, ... % Ignore constraints. 'Fixed', ... {false(4, 1) false(11, 1) false true true true(3, 1) true true});
我们还以同样的方式定义了模型结构的 5 个状态：
InitialStates = struct('Name', {'Angular velocity around x-axis' ... 'Angular velocity around y-axis' ... 'Angular velocity around z-axis' ... 'Angle of attack' 'Angle of sideslip'}, ... 'Unit', {'rad/s' 'rad/s' 'rad/s' 'rad' 'rad'}, ... 'Value', {0 0 0 0 0}, ... 'Minimum', {-Inf -Inf -Inf -Inf -Inf},... % Ignore constraints. 'Maximum', {Inf Inf Inf Inf Inf},... % Ignore constraints. 'Fixed', {true true true true true});
现在，将模型文件与阶数、参数和初始状态数据结合起来，用于创建一个描述该系统的 IDNLGREY 对象：
FileName = 'aero_c'; % File describing the model structure. Order = [5 10 5]; % Model orders [ny nu nx]. Ts = 0; % Time-continuous system. nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ... 'Name', 'Model', 'TimeUnit', 's');
接下来，我们使用 IDDATA 对象中的数据来指定系统的输入和输出信号：
nlgr.InputName = z.InputName; nlgr.InputUnit = z.InputUnit; nlgr.OutputName = z.OutputName; nlgr.OutputUnit = z.OutputUnit;
因此，我们有一个具有 10 个输入信号、5 个状态和 5 个输出信号的 IDNLGREY 对象。如前所述，该模型还包含 23 个参数，其中 7 个为固定参数，16 个为自由参数：
nlgr
nlgr =
Continuous-time nonlinear grey-box model defined by 'aero_c' (MEX-file):
dx/dt = F(t, x(t), u(t), p1, ..., p8)
y(t) = H(t, x(t), u(t), p1, ..., p8) + e(t)
with 10 input(s), 5 state(s), 5 output(s), and 16 free parameter(s) (out of 23).
Name: Model
Status:
Created by direct construction or transformation. Not estimated.
Model Properties
初始模型的性能
在对 16 个自由参数进行估计之前，我们使用初始参数向量对该系统进行仿真。仿真可以提供有关初始模型质量的有用信息：
clf
compare(z, nlgr); % simulate the model and compare the response to measured values
图 3：实际测量结果与初始模型的仿真结果之间的比较。
从图中可以看出，除了 4 到 6 秒这段时间段外，实测信号与仿真信号非常吻合。这一事实在预测误差的绘图中得到了清晰的体现：
figure; h_gcf = gcf; Pos = h_gcf.Position; h_gcf.Position = [Pos(1), Pos(2)-Pos(4)/2, Pos(3), Pos(4)*1.5]; pe(z, nlgr);
图 4：初始模型的预测误差。
参数估计
如上所述，该初始模型是进行参数估计的一个合理的起点。接下来，我们计算这 16 个自由参数的预测误差估计值。此计算需要一些时间。
duration = datetime('now'); nlgr = nlgreyest(z, nlgr, nlgreyestOptions('Display', 'on')); duration = minus(datetime('now'), duration);
估计空气动力学飞行体模型的性能
在所使用的计算机上，参数估计耗时如下：
fprintf('Estimation time : %3.4f seconds\n',seconds(duration))Estimation time : 22.3468 seconds
fprintf('Time per iteration: %3.4f seconds\n',seconds(duration)/nlgr.Report.Termination.Iterations);Time per iteration: 1.0641 seconds
为了评估所估计模型的质量，并说明其相对于初始模型的改进情况，我们对所估计的模型进行仿真，并将实测输出与仿真输出进行比较：
clf compare(z, nlgr);
图 5：实测输出与估计模型的仿真输出之间的比较。
该图清楚地表明，与使用初始模型获得的仿真结果相比，该结果有所改进。现在，对 4 至 6 秒时间段内的系统动态特性的捕捉精度比以前高得多。通过观察预测误差可以最清楚地说明这一点：
figure; h_gcf = gcf; Pos = h_gcf.Position; h_gcf.Position = [Pos(1), Pos(2)-Pos(4)/2, Pos(3), Pos(4)*1.5]; pe(z, nlgr);
图 6：估计模型的预测误差。
让我们通过展示模型和估计的不确定性来结束本案例研究：
present(nlgr);
nlgr =
Continuous-time nonlinear grey-box model defined by 'aero_c' (MEX-file):
dx/dt = F(t, x(t), u(t), p1, ..., p8)
y(t) = H(t, x(t), u(t), p1, ..., p8) + e(t)
with 10 input(s), 5 state(s), 5 output(s), and 16 free parameter(s) (out of 23).
Inputs:
u(1) Aileron angle(t) [rad]
u(2) Elevator angle(t) [rad]
u(3) Rudder angle(t) [rad]
u(4) Dynamic pressure(t) [kg/(m*s^2)]
u(5) Velocity(t) [m/s]
u(6) Measured angular velocity around x-axis(t) [rad/s]
u(7) Measured angular velocity around y-axis(t) [rad/s]
u(8) Measured angular velocity around z-axis(t) [rad/s]
u(9) Measured angle of attack(t) [rad]
u(10) Measured angle of sideslip(t) [rad]
States: Initial value
x(1) Angular velocity around x-axis(t) [rad/s] xinit@exp1 0 (fixed) in [-Inf, Inf]
x(2) Angular velocity around y-axis(t) [rad/s] xinit@exp1 0 (fixed) in [-Inf, Inf]
x(3) Angular velocity around z-axis(t) [rad/s] xinit@exp1 0 (fixed) in [-Inf, Inf]
x(4) Angle of attack(t) [rad] xinit@exp1 0 (fixed) in [-Inf, Inf]
x(5) Angle of sideslip(t) [rad] xinit@exp1 0 (fixed) in [-Inf, Inf]
Outputs:
y(1) V(x)(t) [rad/s]
y(2) V(y)(t) [rad/s]
y(3) V(z)(t) [rad/s]
y(4) Accel.(y)(t) [m/s^2]
y(5) Accel.(z)(t) [m/s^2]
Parameters: Value Standard Deviation
p1(1) Aerodynamic force coefficient [1/(rad*m^2..] 21.2863 0.339394 (estimated) in [-Inf, Inf]
p1(2) -7.62502 0.180264 (estimated) in [-Inf, Inf]
p1(3) 35.0799 0.657227 (estimated) in [-Inf, Inf]
p1(4) 8.58246 1.08611 (estimated) in [-Inf, Inf]
p2(1) Aerodynamic momentum coefficient [1/rad, 1/(..] -1.0476 0.0733533 (estimated) in [-Inf, Inf]
p2(2) 15.6854 0.883102 (estimated) in [-Inf, Inf]
p2(3) 3.00613 0.199227 (estimated) in [-Inf, Inf]
p2(4) -17.7963 0.324639 (estimated) in [-Inf, Inf]
p2(5) -1060.91 224.269 (estimated) in [-Inf, Inf]
p2(6) -53.5594 1.25436 (estimated) in [-Inf, Inf]
p2(7) 34.6095 1.37299 (estimated) in [-Inf, Inf]
p2(8) -210.237 7.95211 (estimated) in [-Inf, Inf]
p2(9) -2641.55 273.034 (estimated) in [-Inf, Inf]
p2(10) -33.6327 3.05742 (estimated) in [-Inf, Inf]
p2(11) -50.9269 1.64086 (estimated) in [-Inf, Inf]
p3 Aerodynamic compensation factor [1/(s*rad)] -0.640669 0.706338 (estimated) in [-Inf, Inf]
p4 Body diameter [m] 0.17 0 (fixed) in [-Inf, Inf]
p5 Body reference area [m^2] 0.0227 0 (fixed) in [-Inf, Inf]
p6(1) Moment of inertia, x-y-z [kg*m^2, kg..] 0.5 0 (fixed) in [-Inf, Inf]
p6(2) 110 0 (fixed) in [-Inf, Inf]
p6(3) 110 0 (fixed) in [-Inf, Inf]
p7 Mass [kg] 107 0 (fixed) in [-Inf, Inf]
p8 Feedback gain [-] 6 0 (fixed) in [-Inf, Inf]
Name: Model
Status:
Termination condition: Maximum number of iterations or number of function evaluations reached..
Number of iterations: 21, Number of function evaluations: 22
Estimated using Solver: ode45; Search: lsqnonlin on time domain data "Data".
Fit to estimation data: [52.93;94.91;91.4;96.07;98.84]%
FPE: 4.627e-10, MSE: 1.672
More information in model's "Report" property.
Model Properties
结束语
该估计模型是研究不同控制策略基本性能的一个良好起点。例如，那些最好具有物理意义的高保真模型，是所谓“模型预测控制系统”的关键组件。
