建模问题学习笔记——经典辨识法

一、数学模型的分类

数学模型是实际物理系统各个物理量之间相互关系的数学描述。

一般数学模型分为线性或非线性、定常或时变、静态或动态、连续或离散、实时或非实时等形式。

1、线性或非线性的区别

$$a_0\frac{d^n}{d{t}^n}c(t)+a_1\frac{d^{n-1}}{d{t}^{n-1}}c(t)+...+a_{n-1}\frac{d}{dt}c(t)+a_{n}c(t)=b_0\frac{d^m}{d{t}^m}r(t)+b_1\frac{d^{m-1}}{d{t}^{m-1}}r(t)+...+b_{m-1}\frac{d}{dt}r(t)+b_{m}r(t)$$
可用以上形式的微分方程描述的系统则为线性系统。
$$\ddot{y}(t)+y(t)\dot{y}(t)+y^2(t)=r(t)$$
类似以上形式的微分方程描述的系统则为非线性系统，其最主要的特征是含有变量及其导数的高次幂或乘积项或者说是系数与变量有关。

2、静态或动态的区别

一般定常或时变描述线性系统，因此针对线性系统的微分方程描述：$$a_0\frac{d^n}{d{t}^n}c(t)+a_1\frac{d^{n-1}}{d{t}^{n-1}}c(t)+...+a_{n-1}\frac{d}{dt}c(t)+a_{n}c(t)=b_0\frac{d^m}{d{t}^m}r(t)+b_1\frac{d^{m-1}}{d{t}^{m-1}}r(t)+...+b_{m-1}\frac{d}{dt}r(t)+b_{m}r(t)$$
若其系数是常数，则该线性系统为线性定常系统；若其系数随时间而变化，则该线性系统为线性时变系统。

3、连续或离散的区别

离散系统的信号在时间上是离散的，因此一般用差分方程来描述。连续信号通过采样开关的采样就可以转换为离散信号。

一般离散系统的表现形式为：
$$a_{0}c(k+n)+a_{1}c(k+n-1)+...+a_{n-1}c(k+1)+a_{n}c(k)=b_{0}r(k+m)+b_{1}r(k+m-1)+...+b_{m-1}r(k+1)+b_{m}r(k)$$

二、建立模型的方法

一般建立数学模型的方法有解析法和试验法。

解析法是通过其内部运动规律等，利用有关定理，用数学方法建立数学模型。

试验法是通过对输入输出数据进行处理，获取其特性从而得到数学模型的方法，也称为系统辨识。

三、试验法建模

系统辨识方法建立起的被控对象的数学模型可分为参数模型（如传递函数、状态空间方程等）和非参数模型（如阶跃响应、脉冲响应），非参数模型可以通过变换转化为参数模型。

系统辨识的方法可以分为经典和现代两大类。经典辨识法有时域法、频率响应法、相关辨识法等，现代辨识法有最小二乘法、极大似然法、随机逼近法等。

1、时域法

时域法辨识对象的动态特性是用试验方法测出对象的响应曲线，与标准的传递函数响应曲线进行比较，确定所辨识的对象属于哪一类传递函数，并从响应曲线中求出传递函数的各个参数。常用的方法有阶跃响应法和脉冲响应法。

1.1 阶跃响应法

阶跃响应法是当系统处于稳态时，快速施加一个阶跃扰动，并记录其输入输出数据。

1.2 脉冲响应法

脉冲响应法是当系统处于稳态时，快速施加一个脉冲扰动，一般用于某系统无法适用于长时间阶跃输入时。脉冲施加方法为对被控对象施加一个阶跃扰动，直到被控对象即将到达允许的最大偏差值时再清除扰动，脉冲宽度可以进行适当的选择。

测定的脉冲响应曲线需要换算成阶跃响应曲线才能与标准传递函数的响应曲线进行比较。

(1)脉冲响应曲线到阶跃响应曲线的换算

由于矩形脉冲可以等效为两个阶跃信号的叠加，因此脉冲响应曲线可以看成是0时刻开始的一条阶跃响应曲线和t0时刻开始的一条阶跃响应曲线的叠加：

1. 我们使用脉冲响应法得到了OADE这样一条脉冲响应曲线；
2. 由于上下两个叠加的阶跃信号大小相等，方向相反。所以FG段曲线应当和OA段曲线相同，我们可以通过AD段曲线和FG段曲线得到AB段曲线；
3. 同理可以GH段曲线和AB段曲线是相同的，得到GH段曲线和DE段曲线得到BC段曲线。之后亦同理。

1.3 阶跃响应曲线求传递函数

(1)一阶惯性环节系统辨识

一阶惯性环节系统的阶跃响应曲线如上所示。
一阶惯性环节系统的传递函数描述为：$$G(s)=\frac{K}{Ts+1}$$
一阶惯性环节系统的阶跃响应曲线数学描述为：$$y^{\ast }(t)=1-e^{-t/T}$$
当$t=T$时，$y^{\ast }(t)=0.63$；当$t=2T$时，$y^{\ast }(t)=0.87$。

施加$U_0$阶跃扰动时，响应曲线幅值为$y(\infty )$，则系统的增益为$K=y(\infty )/U_0$，并通过响应曲线的数学描述求取时间常数$T$。

(2)带纯时延的一阶惯性环节——S形曲线

带纯时延的一阶惯性环节系统的阶跃响应曲线如上所示。
带纯时延的一阶惯性环节系统的传递函数描述为：$$G(s)=\frac{K}{Ts+1}{e}^{-\tau s}$$
一阶惯性环节系统的阶跃响应曲线数学描述为：$$y^{\ast }(t)=\{\begin{array}{ll}0,& t<\tau \\ 1-e^{-(t-\tau )/T},& t\ge \tau \end{array}$$
当$t=t_1$时，$y^{\ast }(t)=0.39$；当$t=t_2$时，$y^{\ast }(t)=0.63$。而$T=2(t_2-t_1)$，$\tau =2t_1-t_2$。

施加$U_0$阶跃扰动时，响应曲线幅值为$y(\infty )$，则系统的增益为$K=y(\infty )/U_0$，并通过响应曲线的数学描述求取时间常数$T$和时延$\tau$。

2、频率法

通过在系统的输入端施加不同频率的正弦信号，然后测量系统输出的稳态响应，再根据幅值比和相位差作出系统的对数频率特性曲线。再将对数幅频曲线用斜率为0dB/dec、$\pm 20$dB/dec、$\pm 40$dB/dec…直线分段近似，获得对数幅频渐进特性曲线，并确定最小相位条件下系统的传递函数。

1. 确定系统积分或微分环节的个数。对数幅频渐进特性曲线的低频渐近线的斜率为$-20\nu$dB/dec，$\nu$正为积分环节，否则为微分环节。
2. 确定系统传递函数结构形式。根据交接频率斜率变化确定环节的种类，-20dB/dec对应惯性环节，-40dB/dec对应重惯性环节（谐振为振荡环节）。
3. 确定传递函数参数。由低频渐进特性线方程得增益K、交接频率。

3、相关辨识法

白噪声信号对系统的正常运行状态影响不大，不会使被测对象偏离正常工作状态，但是白噪声的产生比较困难，因此常采用伪随机信号作为相关辨识法的输入信号，即M序列。

3.1 M序列的产生

M序列是用n个移位寄存器所能得到的周期最长的二元序列信号，是伪随机二元序列信号的一种。

移位寄存器的个数为n，M序列一个周期的拍数为$N=2^n-1$，模二门就是两个移位寄存器的输出作异或运算（相同为0，不同为1）。

%% 产生M序列
r=0.9;
a=1.1;
clc;clear all
N=31;
Amplitude=0.5;offset=3;
L=10*N;% M序列的周期
ym1=1;ym2=1;ym3=1;ym4=0;ym5=0;%五个移位积存器的输出初始值
for i=1:L
    ym1_=xor(ym5,ym3);      % 异或 相同为0 不同为1
    ym2_=ym1;
    ym3_=ym2;
    ym4_=ym3;
    ym5_=ym4;
    ym(i)=ym5;
    if ym(i)>0.5
        u(i)=-Amplitude+offset;            % M序列的幅值为5
    else
        u(i)=Amplitude+offset;
    end
    ym1=ym1_;ym2=ym2_;ym3=ym3_;ym4=ym4_;ym5=ym5_;
end 
% 第一种含噪声模型
v = randn(1,L)*0.001;
y(1) = 0+v(1);t(1)=0;
for i=2:1:L
    t(i)=t(i-1)+0.01;
    y(i) = 0.8187*y(i-1)+0.1876*u(i-1)+v(i);
end
t=t';u=u';y=y';
M_list=[t u y];
save('M_list.mat','M_list');

3.2 用M序列辨识系统

以M序列作为输入信号，输入到系统得到输出信号，然后通过输入输出信号求它们的互相关函数，进一步可以求出系统的脉冲响应函数。

clear
clc
delta=0.01;N=31;
a=0.5;  %输入信号x的幅值
%% 导入M序列的输入输出数据并进行处理
load('M_list.mat');
t=M_list(100:300,1)';%取符合周期性规律的数据
X=M_list(100:300,2)';
Y=M_list(100:300,3)';
X=X-3;  % 减去偏置量
y0=sum(Y)/(300-100+1);
Y=Y-y0;
t=t-1;
%% 计算互相关函数（一个周期）
Rxy=zeros(N,1);
for i=0:(N-1)
    sum=0;
    for k=0:(N-1)
        xid=interp1(t,X,k*delta);
        yidtao=interp1(t,Y,k*delta+i*delta);
        sum=sum+xid*yidtao;
    end
    Rxy(i+1)=sum/N;
end
%% 根据互相关函数求脉冲响应
Rxy=Rxy-min(Rxy);
gt=N/(a^2)/(N+1)*Rxy/delta;
figure(1);
plot(0:0.01:(N-1)*0.01,gt,'r');
%% 根据脉冲响应绘制阶跃响应曲线，从而辨识系统 
stepResponse=zeros(N,1);
for i=1:N
    for j=1:i
        stepResponse(i)=stepResponse(i)+delta*gt(j);
    end
end
figure(2);
plot(0:0.01:(N-1)*0.01,stepResponse,'b');
%% 辨识系统
K=(stepResponse(N)-stepResponse(1));
y_T=0.63*(stepResponse(N)-stepResponse(1));
for i=1:1:50
    if stepResponse(i+1)>y_T
        T=(y_T-stepResponse(i))/(stepResponse(i+1)-stepResponse(i))*delta+(i-1)*delta;
        break;
    end
end
model_es=tf([K],[T,1])

本文是作者在日常学习生活中所作，难免有遗漏或错误，遇到问题的读者请点击给我写信向我的邮箱反馈，不胜感激。