MATLAB 非线性系统辨识（电子节气门系统辨识）

MATLAB 非线性系统辨识（电子节气门系统辨识）

最近筹备新的文章，主题是节气门的固定时间滑模控制。在实验的时候系统标称参数和实际系统相去甚远，做跟踪控制时总有一段稳态误差。遂重新辨识一下系统的参数。

此前从未辨识过系统参数，参考了 MATLAB 帮助文档中的例子，对节气门系统进行了辨识。

节气门（ET）系统建模

考虑节气门体的弹簧非线性、摩擦非线性可以得到一个简洁的系统模型：
$$J_{eq}{\ddot{\theta }}_t+B_{eq}{\dot{\theta }}_t+F_s\mathrm{sign}{\dot{\theta }}_t+{\tau }_{LH}\mathrm{sign}({\theta }_t-{\theta }_0)+k_{sp}({\theta }_t-{\theta }_0)+T_d=bu$$

简单整理一下，可以得到回归模型：

$${\ddot{\theta }}_t=a_{1}u-a_2\mathrm{sign}{\dot{\theta }}_t-a_3{\dot{\theta }}_t-a_4({\theta }_t-{\theta }_0)-a_5\mathrm{sign}({\theta }_t-{\theta }_0)-w$$

输入输出数据收集

节气门电机的额定电压是 $12V$，这里为了不让节气门过快的饱和，使用了 0 到 3V的正弦波信号，周期 $T$ 是 $1s$ 和 $0.5s$：
$$u=1.5+1.5\sin (2\pi /T\cdot t)$$
采集输入信号 $u$，以及输出信号 ${\theta }_t$（阀片的转动弧度位置）用于辨识。

MATLAB nlgrey 函数辨识非线性系统

简单列出步骤：

定义 ODE 函数

令 $x_1={\theta }_t,x_2={\dot{\theta }}_t,w=0$, 将回归模型改写为状态方程形式：
$$\begin{gathered} \{\begin{array}{rl}{\dot{x}}_1& =x_2\\ {\dot{x}}_2& =a_{1}u-a_2\mathrm{sign}{x}_2-a_3{x}_2-a_4(x_1-{\theta }_0)-a_5\mathrm{sign}(x_1-{\theta }_0)\end{array} \\ y=x_1 \end{gathered}$$
新建 ET_m.m 文件，定义节气门 ODE 函数(m 文件形式)用于辨识

function [dx,y] = ET_m(t,x,u,a1,a2,a3,a4,a5,varargin)
th0 = 0;
y    = x(1);
dx  = [x(2)
   	     a1 * u - a2 * sign(x(2)) - a3 * x(2) - a4 * (x(1) - th0) - a5 * sign(x(1) - th0)];
end

定义主函数

新建 ET_identification.m 文件，进行辨识

idnlgrey 函数定义一个用于回归的灰箱模型：

Note:

• 参数的初值参考自己需要辨识的系统的相关论文
• 起名字可以随意，不必按照本例

%--- System model def ---

clc,clear,close all
FileName      = 'ET_m';
Order         = [1 1 2];     % Order of Output, input and state, [ny nu nx];
th0           = 0;
Parameters    = [1400;94;68;28;186];  % Parameters initilization
InitialStates = [th0;0];    % Initial state
Ts            = 0; % Sampling time, 0 means continous system
ETsys = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ...
                'Name', 'ET Model', ...
                'InputName', {'Voltage'}, ...
                'InputUnit', {'V'}, ...
                'OutputName', {'Position'},...
                'OutputUnit', {'rad'}, ...
                'TimeUnit', 's');

定义系统状态

%--- System state def ---

ETsys = setinit(ETsys, 'Name', {'Position';'Velocity'});
ETsys = setinit(ETsys, 'Unit', {'rad'; 'rad/s'});

定义待辨识参数

• 如果某些参数已知，可以将其固定，参数的范围最好限定一下，本例没有做限定

%--- Identified name def ---
ETsys = setpar(ETsys, 'Name', {'a1'; 'a2'; 'a3'; 'a4';'a5'}); % give a readible name for the parameters
% ETsys.Parameters(1).Fixed = true; % If some parameters are fixed, pls declare here
ETsys.Parameters(1).Minimum = 0;
ETsys.Parameters(1).Maximum = inf;
ETsys.Parameters(2).Minimum = 0;
ETsys.Parameters(2).Maximum = inf;
ETsys.Parameters(3).Minimum = 0;
ETsys.Parameters(3).Maximum = inf;
ETsys.Parameters(4).Minimum = 0;
ETsys.Parameters(4).Maximum = inf;
ETsys.Parameters(5).Minimum = 0;
ETsys.Parameters(5).Maximum = inf;

present(ETsys)  % Print system information

导入开环响应数据，建立实验数据模型

• 本例中输入量归一化不是必要的，根据自己的系统来

%---- Input experimental data ---

ETdata         = csvread('ModelIdentification3.csv',3,0); % Input data, my data start from row 3, column 0
ns             = 1666;   % start point, my data start from 16s
nf             = 2300; % end point, my data end in 23s
dataset.Input  = ETdata(ns:nf,2) ./ 12;      % u, normalized input value
dataset.Output = ETdata(ns:nf,4); % thetat is in row 4 in my case
th0            = 0;       % LH position for ET system
SampleTime     = .01; % Sampling time

z            = iddata(dataset.Output, dataset.Input, SampleTime, 'Name', 'Actual Data');
z.InputName  = ETsys.InputName;
z.InputUnit  = ETsys.InputUnit;
z.OutputName = ETsys.OutputName;
z.OutputUnit = ETsys.OutputUnit;
z.Tstart     = 0;
z.TimeUnit   = 's';
present(z); % print print experimental system info

开始辨识

最主要的是 nlgreyest 函数，其他都是展示结果的语句

X0init = [th0;0]; % Initial states
ETsys  = setinit(ETsys, 'Value', num2cell(X0init)); 
figure
compare(z, ETsys, compareOptions('InitialCondition', X0init)); % comparsion of experimental data and the model with initial parameters
ETsys  = nlgreyest(ETsys, z, nlgreyestOptions('Display', 'on')); % Starting identification
figure
compare(z, ETsys, compareOptions('InitialCondition', X0init)); % comparson of experimental data with the model under identified parameters
present(ETsys) % print identified system info

辨识效果

初始参数下的数据和模型对比：

辨识过程，这个过程比较长，本例用了 5 mins 左右，cost 至少小于 $0.001$ 才算差不多：

辨识后的系统和实验数据对比：

辨识出来的参数：

使用新的参数之后，跟踪效果大大提高。