智能机器人与运动控制-第11讲-预报误差极小化辨识的Matlab实现

智能机器人与运动控制11.预报误差极小化辨识的Matlab实现研究背景与基础概念回顾MatlabIdgrey方法与灰箱建模系统辨识程序主体仿真数据生成与验证研究背景与基础概念回顾研究背景与意义直流电机的工业应用：驱动系统、精密控制等场景。模型辨识的重要性：为控制算法设计提供准确数学模型预报误差极小化：通过最小化“实际输出与模型预报输出的差异”，实现参数最优估计。核心内容回顾

递推最小二乘估计法

MatlabIdgrey方法与灰箱建模MatlabIdgrey方法介绍

状态空间灰箱模型构建——函数定义function[A,B,C,D,K1]=motor_model(par,type,T)J=1.0258e-06;Ce=0.002593;R=3.4453779;A=[010;-par(1)/J-par(2)/Jpar(3)/J;0-Ce*60*par(4)/(2*pi)-R*par(4)];B=[0;0;par(4)];

代码%输入输出：par（待辨识参数）、type（模型类型）、T（采样时间）；输出A,B,C,D（状态空间矩阵）与K1（Kalman滤波器）。%依据之前拟合得到%连续性形式的含参A矩阵（物理对应：par(1)=K（转矩系数）、par(2)=Bv（粘滞摩擦系数）、par(3)=Cm（电磁转矩系数）、par(4)=1/L（电感倒数））%连续性形式的含参B矩阵解释模型离散化处理AB=expm([AB;zeros(1,4)]*T);

A=AB(1:3,1:3);

B=AB(1:3,4);C=[0,1,0;0,0,1];D=[0;0];

代码%依据采样间隔离散化模型%提取离散A矩阵%提取离散B矩阵\%C矩阵%D矩阵解释模型离散化处理R1=zeros(3);R2=1e-4*eye(2);[~,K1]=kalman(ss(A,eye(3),C,zeros(2,3),T),R1,R2);

代码%过程噪声、观测噪声%将状态空间模型转换为“输入-输出”形式，便于辨识算法处理解释系统辨识程序主体辨识程序主体——参数初始化theta_1=[0.001;0.01;257.4;0.01696];

dt=0.01;

代码%待辨识参数初始值%采样时间解释辨识程序主体——idgrey模型生成Minit=idgrey('motor_model',theta_1,'d',dt);Minit.Structure.Parameters.Minimum=[-Inf;Inf;Inf;Inf];

Minit.Structure.Parameters.Maximum=[Inf;Inf;Inf;Inf];data_id=iddata(y,u,dt);

代码%生成Matlab灰箱模型%参数下界%参数上界%封装输出y、输入u、采样时间dt，将原始数据转换为Matlab辨识工具箱可识别的iddata对象解释辨识程序主体——idgrey模型生成opt=greyestOptions;

opt.InitialState='estimate';

opt.Display='full';

opt.SearchOptions.MaxIterations=20;model_est=greyest(data_id,Minit,opt);

代码%设置辨识参数%初始状态设为“待估计”%显示完整辨识过程%最大迭代次数%执行灰箱模型辨识解释仿真数据生成与验证仿真数据生成——目的与意义为什么需要仿真数据？缺乏真实实验数据时，验证辨识算法的正确性；模拟“带噪声的真实场景”，测试算法抗干扰能力。仿真数据生成——真实参数设置K=0.06;J=1.0258e-06;Bv=4.27;Cm=0.32/5.9;Ce=0.002593;L=1.12e-3;R=3.453779;dt=0.01;

代码%代码（直流电机真实物理参数）解释仿真数据生成——状态空间矩阵（真值）A=[010;-K/J-Bv/JCm/J;0-Ce*60/(L*2*pi)-R/L];B=[0;0;1/L];C=[1,0,0;0,0,1];D=zeros(2,1);

代码%A矩阵%B矩阵%C矩阵%D矩阵解释仿真数据生成——输入与噪声序列n=50000;u=50*randn(n,1);v=sqrt(1e-5)*randn(1,n);

代码%输入序列（模拟激励）：%数据个数%随机输入（正态分布）%噪声序列（模拟测量噪声）：%随机噪声解释仿真数据生成——状态转移与输出计算x_1=zeros(3,1);

AB=expm([AB;zeros(1,4)]*dt);A=AB(1:3,1:3);B=AB(1:3,4);y=zeros(2,n);x=zeros(3,n);fori=1:n-1x(:,i+1)=A*x(:,i)+B*u(i);y(:,i)=C*x(:,i)+v(:,i);endy(:,n)=C*x(:,n)+v(:,n);x=x';y=y'代码%初始状态%初始化输出、状态%状态转移%输出计算（加噪声）%最后一个输出%转置为行向量，适配辨识数据格式解释辨识结果展示——输出误差直方图误差集中在0附近，呈正态分布，说明辨识模型