利用辨识工具箱辨识状态空间表达式 | 学习切片

最新推荐文章于 2024-06-07 11:16:18 发布

镭射麒麟臂

最新推荐文章于 2024-06-07 11:16:18 发布

阅读量1k

点赞数 21

文章标签： matlab 开发语言

本文链接：https://blog.csdn.net/weixin_52059170/article/details/134647058

版权

本文介绍了MATLAB中用于线性（idgrey）和非线性（idnlgrey）模型辨识的函数，包括如何处理状态空间表达式的参数估计，提供线性模型的示例和非线性模型在实际问题中的应用实例，以及如何设定参数约束以提高辨识精度。

摘要由CSDN通过智能技术生成

灰盒估计

对于状态空间表达式 $\left\{\begin{matrix} \dot{x}=Ax+Bu\\y=Cx+Du \end{matrix}\right.$ ，其中，u是输入，y是输出，x是状态变量，A，B，C，D均未知，其中A为 $N\times N$ 矩阵，B为 $N\times 1$ 矩阵，C为 $N\times N$ 矩阵，D为 $N\times 1$ 矩阵。

与黑盒估计不同，这里的A，B，C，D矩阵之间包含约束关系，如果使用一般的辨识函数如ssest或n4sid，无法保证辨识过程中始终保持约束关系，因此我们这里需要使用其他函数。

下面我们针对线性模型和非线性模型的辨识采用的辨识函数的使用进行介绍。

一、线性模型的辨识

对于线性模型 $\left\{\begin{matrix} \dot{x}=Ax+Bu\\y=Cx+Du \end{matrix}\right.$ ，我们使用到名为idgrey的函数，即为grey-box model with identifiable parameters，表示灰盒模型的线性参数辨识。

函数功能

idgrey函数在matlab中使用的基本形式为：

sys = idgrey(odefun,parameters,fcn_type,extra_args,Ts)

其中，输入的参数表示的含义如下：

输入参数	含义
odefun	需要辨识的ODE函数或MATLAB中已有的ODE函数模型
parameters	对需要辨识的参数的初始值进行定义
fcn_type	模型是在连续时间、离散时间还是同时参数化； 'c'表示返回连续时间对应矩阵，与Ts无关； 'd'表示返回离散时间对应矩阵，与Ts可能有关也可能无关； 'cd'表示如果Ts=0，为连续时间；Ts>0，为离散时间
extra_args	odefun所需输入的额外参数（如没有可以不写）
Ts	采样时间连续时间，Ts=0；离散时间，Ts=-1；具有指定采样时间的离散模型，用指定的TimeUnit表示

程序实现

示例一

下面以一个简单的例子演示模型的使用方法。

对于表达式 $\dot{x}(t)=\begin{bmatrix} 0 & 1\\ 1 & -\frac{1}{\tau} \end{bmatrix}x(t)+\begin{bmatrix} 0\\ \frac{G}{\tau} \end{bmatrix}u(t)$ ， $y=\begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix}x(t)$ ，其中 $G$ 为已知参数， $\tau$ 为要估计的参数，方程满足状态空间表达式 $\left\{\begin{matrix} \dot{x}(t)=Ax(t)+Bu(t)\\y(t)=Cx(t) \end{matrix}\right.$ ，可以在matlab中进行估计。

第一步，将表达式写为ODE函数：

function [A,B,C,D] = fcn(tau,ts,G)
A = [0 1;0 -1/tau];
B = [0;G/tau];
C = eye(2);
D = [0;0];
end

这里，tau是要辨识的参数；G为已知参数，数值在extra_args中给出即可，也可以直接输入；ts为采样时间。

第二步，给出idgrey函数的输入：

odefun = 'fcn';
parameters = 1;
fcn_type = 'cd';
extra_args = 0.25;
Ts = 0;

fcn就是第一步给出的函数；parameters是要辨识的参数tau的初始值=1；选择的函数类型为'cd'，且Ts=0，说明这里采用的是连续时间；extra_args=0.25，其实就是G=0.25。

第三步，利用idgrey函数进行辨识：

sys = idgrey(odefun,parameters,fcn_type,extra_args,Ts)

辨识后可以得到参数tau=1，在工作区点击sys可以查看A,B,C,D的值。

示例二

对于有更多约束的模型，可以添加一些限定条件，如固定一些参数或者对参数变化范围进行限定，这样可以提高辨识的准确性和效率，也可以避免在实际应用中，迭代出现局部收敛的情况。

下面举一个简单的示例来体现这种操作。

对于表达式 $\dot{x}=\begin{bmatrix} 0 & 1\\ -\frac{g}{l} & -\frac{b}{ml^{2}} \end{bmatrix}x$ ， $y=\begin{bmatrix} 1&0 \end{bmatrix}x$ ，其中，m为模型质量，g为重力加速度，l为模型长度，b为粘性摩擦系数，下面对该模型的参数进行辨识。

第一步，将表达式写为ODE函数：

function [A,B,C,D] = LinearPendulum(m,g,l,b,Ts)
A = [0 1; -g/l, -b/m/l^2];
B = zeros(2,0);
C = [1 0];
D = zeros(1,0);
end

第二步，给出idgrey函数的输入：

odefun = 'LinearPendulum';
m = 1;
g = 9.81;
l = 1;
b = 0.2;
parameters = {'mass',m;'gravity',g;'length',l;'friction',b};
fcn_type = 'c';

LinearPendulum就是第一步给出的函数；parameters是要辨识的参数m,g,l,b的初始值；选择的函数类型为'c'，说明这里采用的是连续时间。

第三步，针对辨识参数的物理意义，对它们进行约束：

sys.Structure.Parameters(1).Free = false;
sys.Structure.Parameters(2).Free = false;
sys.Structure.Parameters(3).Free = false; 
% sys.Structure.Parameters(2).Fixed = 1;%0;
sys.Structure.Parameters(4).Minimum = 0; 
% sys.Structure.Parameters(4).Maximum = 0;

这里表示，第1,2,3个参数在迭代过程中固定不变，第4个参数在迭代过程中的最小值为0，即：m,g,l始终不变，b>0。

如果想表示参数在迭代过程中变化，可以写为Free = true或Fixed = 0；Fixed = 1表示参数自由变化。Maximum和Minimum表示对参数最大值和最小值的约束。

第四步，利用idgrey函数进行辨识：

sys = idgrey(odefun,parameters,fcn_type);

辨识结果可以在parameters中查看。

二、非线性模型的辨识

当A,B,C,D中的参数为时变时，上述的微分方程就变为了非线性模型，如果仍采用上述函数进行辨识，会产生很大的误差，因此我们需要使用MATLAB中的非线性辨识函数idnlgrey来辨识。

函数功能

idnlgrey函数在matlab中使用的基本形式为：

sys = idnlgrey(FileName,Order,Parameters,InitialStates,Ts)

其中，输入的参数表示的含义如下：

输入参数	含义
FileName	存储模型结构的函数名称
Order	模型的输出、输入和状态数 [Ny,Nu,Nx]表示模型输出的数量Ny，输入的数量Nu，状态数Nx
Parameters	模型参数，可以包含以下信息：名称，单位，初值；也可以对参数进行约束
InitialStates	辨识参数的初始值
Ts	采样时间

程序实现

下面用一个具体示例展示idnlgrey函数的使用方法，使用的是MATLAB辨识工具箱文档中的一个例子。

对于 $\left\{\begin{matrix} \dot{x}=Ax+Bu\\y=Cx+Du \end{matrix}\right.$ ，其中，A,B,C,D矩阵中含有时变的参数，我们使用辨识方法得到模型的动力系数和力矩系数。

第一步，我们导入测量数据：

load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'aerodata'));

使用的是MATLAB工具箱现成的数据，在自己使用时导入自己测量的数据。

第二步，将测量数据利用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';

这里，步长dt为0.02s，u是输入，y是输出；输入u包含：俯仰舵偏角，滚转舵偏角，偏航舵偏角，动压，速度，x-y-z三轴的角速度，攻角，侧滑角；输出y包含：x-y-z三轴速度，y-z两轴加速度。

第三步，得到输入数据的图像：

figure('Name', [z.Name ': input data'],...
   'DefaultAxesTitleFontSizeMultiplier',1,...
   'DefaultAxesTitleFontWeight','normal',...
   'Position',[50, 50, 850, 620]);
for i = 1:z.Nu
   subplot(z.Nu/2, 2, i);
   plot(z.SamplingInstants, z.InputData(:,i));
   title(['Input #' num2str(i) ': ' z.InputName{i}],'FontWeight','normal');
   xlabel('');
   axis tight;
   if (i > z.Nu-2)
       xlabel([z.Domain ' (' z.TimeUnit ')']);
   end
end

第四步，得到输出数据的图像：

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:z.Ny
   subplot(z.Ny, 1, i);
   plot(z.SamplingInstants, z.OutputData(:,i));
   title(['Output #' num2str(i) ': ' z.OutputName{i}]);
   xlabel('');
   axis tight;
end
xlabel([z.Domain ' (' z.TimeUnit ')']);

第五步，动力学建模，这里使用的areo.c是空气动力学中的模型，使用状态空间表达式，输出方程表示，MATLAB的help文档可以找到具体描述。

第六步，设置各参数的初值：

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});
% 4个力系数，11个力矩系数，空气动力学补偿系数，翼展，机翼参考面积，转动惯量，质量，反馈增益

在这里可以改变4个力系数和11个力矩系数的初值，以及翼展，参考面积，质量转动惯量这些系数。

第七步，如果需要，可以对力系数和力矩系数的范围进行限制，里面的范围可以自行更改：

Parameters(1).Minimum =[-Inf;-Inf;-Inf;-Inf];
Parameters(1).Maximum =[Inf;Inf;Inf;Inf];
Parameters(2).Fixed = [0;1;0;0;1;0;0;1;1;1;0];
Parameters(2).Maximum = [Inf;Inf;Inf;Inf;Inf;Inf;Inf;Inf;Inf;Inf;Inf];
Parameters(2).Minimum = [-Inf;-Inf;-Inf;-Inf;-Inf;-Inf;-Inf;-Inf;-Inf;-Inf;-Inf];

第八步，定义模型结构的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});

第九步，开始辨识：

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

第十步：将已知的输出向量与模型参数初始值得到的输出向量进行比较：

nlgr.InputName = z.InputName;
nlgr.InputUnit = z.InputUnit;
nlgr.OutputName = z.OutputName;
nlgr.OutputUnit = z.OutputUnit;

nlgr

clf
compare(z, nlgr);

第十一步，评估估计参数的准确性，将利用辨识得到的参数计算的输出与真实测量得到的输出数据进行比对：

duration = clock;
nlgr = nlgreyest(z, nlgr, nlgreyestOptions('Display', 'on')); 
duration = etime(clock, duration);

clf
compare(z, nlgr);

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

第十二步，输出结果：