201404703028系统辨识资料

1、玄倨科技天号研究生课程实验封面课程名称：报告题目：学生学号：学生姓名：任课教师：学位类别：评分标准 及分值选题与参阅资料（分值10）报告内容（分值 60 ）报告表述（分值20 ）创新性（分值 10 ）评分报告评语：总评分评阅教师：评阅时间年 月曰注：此表为每个学生的报告封面，请任课教师填写分项分值1系统辨识部分1.1 题目直流电动机的辨识与自适应系统设计与仿真1.1.1 工作原理图1是一台最简单的直流电动机的模型，N和S是一对固定的磁极（一般是 电磁铁，也可以是永久磁铁）。磁极之间有一个可以转动的铁质圆柱体，称为电 枢铁芯。铁芯表面固定一个用绝缘导体构成的电枢线圈abcd,线圈的两端分别接到相

2、互绝缘的两个弧形铜片上，弧形铜片称为换向片，它们的组合体称为换向 器。在换向器上放置固定不动而与换向片滑动接触的电刷A和B,线圈abcd通过换向器和电刷接通外电路。电刷铁芯、电刷线圈和换向器构成的整体称为电枢。此模型作为电动机运行时，将直流电源加于电刷A和B,例如将电源正极加 于电刷A,将电源负极加于电刷B,则线圈abcd中流过电流。在导体ab中，电 流由a流向b,在导体cd中，电流由c流向do载流导体ab和cd均处于N和S 极之间的磁场中，受到电磁力的作用。电磁力的方向由左手定则确定， 可知这一 对电磁力形成一个转矩，称为电磁转矩，转矩的方向为逆时针方向，使整个电刷 逆时针方向旋转。当电刷旋

3、转180度，导体cd和ab交换位置，如图1所示。由 于电流仍从电刷A流入，使cd中的电流变为由d流向c,而ab中的电流由b流 向a,从电刷B流出，用左手定则判别可知，电磁转矩的方向仍是逆时针方向。由此可见，加于直流电动机的直流电流，借助换向器和电刷的作用，变为电枢线圈中的交变电流。这种将直流电流变为交变电流的过程称为逆变。由于电刷线圈所处的磁极也是同时交变的， 这使电刷产生的电磁转矩的方向恒定不变， 从 而确保直流电动机朝确定的方向连续旋转。这就是直流电动机的基本工作原理。图1直流电机工作的基本工作原理1.2 辨识目的对于一个系统，为什么要对它进行辨识？辨识的意义又在哪里？这些问题在 很久以前

4、就有人提出过，也进行过深入的研究，并总结出一套成熟的方法：最小 二乘辨识方法、最大似然辨识方法、梯度法辨识等等。然而，这些方法在线性问 题上容易解决，而在非线性问题上却相对较复杂， 且方法并非唯一，而且找不到 统一的设计模式，只能是针对具体问题分析其分线性的问题所在， 抓住其影响系 统动、静态品质的要害，研究辨识非线性系统模型及控制的理论和方法，进而对 系统进行辨识、补偿或控制。所以，若能够通过辨识得到其精确的模型，则是控制问题的关键。而本设计 中，涉及到的也是相同的问题，对于直流电机，我们给定输入，由输出与输入的 相应关系，我们便可以通过辨识得到直流电机的模型参数， 进行可以进行模型精 确的

5、控制，而避免了对模型的不确定控制。1.3 直流电机的数学模型直流电机电枢回路的电路方程是：di.ua - E = iRa La (1.1)式(1.1)中，Ua是加到电机两端的电压；E是电机反电势；i是电枢电流；Ra是电枢回路总电阻；La是电枢回路总电感；I aTi =上 称为电枢回路电磁时间常数。Ra并且反电动势E与电机角速度6 m成正比：E=K % =/m(1.2)式(1.2)中，ke称为反电势系数；与为电机轴的转角。对于电机而言，其转动轴上的力矩方程为：kmi -M1 = Jmlm =Jm4(1.3)式(1.3)中，km是电机的力矩系数;Ml是负载力矩;Jm是电机电枢的转动惯量。对式（1.

6、1）、（1.2）、（1.3）进行拉氏变换得到:Ua(s)-E(s)=Ra(I(s) +TiI(s)s)E(s) =kJm(s)(1.4)、kmI (s)-Mi =JmL(s)s由式（1.4）,可得从电枢电压Ua到转速4的传递函数:1/ Ra t s _ Tis 1 Ua s 1 . 1/ RaT1s 1J mskmkeJmskm/ Ra2TJms Jms 。/ Ra(1.5)很明显的看出，这是个典型的二阶系统的传递函数1.4 系统的输入与输出由式（1.4）的方程组可以得到相应的电动机数学模型的结构框图:根据图2所示的直流电动机数学模型结构图，我们定义系统的输入与输出分 别为：系统的输入：给定电

7、压Ua系统的输出：电机转速41.5 系统的数据采集图3数据采集直流电动机的卒&入给定电压Ua通过A/D转换给于Mcu,同时直流电动机的 输出转速&经编码器测量得到，将脉冲数送于 Mcu, Mcu再经用口将数据发送 给电脑。此时，电动机系统的输入与输出也得知， 便可通过电脑分析求解出电机 的系统模型。其中，A/D转换采用16位精度，输入电压范围711V, 100Hz采 样频率的芯片；编码器采用500线精度，表示电机每转一圈，输出 500个脉冲。1.6 输入信号的选择输入采用幅值为7V的方波信号，如图4所示。图4方波信号1.7 系统的仿真由于前面所述的传递函数是连续函数，所以必须将

8、其离散化，取采样时间为0.9sos s式(1.4)为包Ua Skm/RaTlJms2 Jms kmke/Ra,现取电机的传递函数为:苗(S)Ua(S)176.20.01S20.1s 1(1.6)对上式做z变换，得到离散化的传递函数为:(1.7)N(z) _59.96 z 42.59U(z) - z2 - 0.7859 z 0.3679将其转换为差分方程为：N(K) =0.7859N(K -1) -0.3679N (K -2) 59.96U (K -1) 42.59U(K -2)当然，一个系统并不都是理想的系统，还包括各种系统干扰，以e(K)表示方差为零的系统干扰，于是差分方程就变为：N(K)

9、=0.7859N(K -1) -0.3679N(K -2) 59.96U (K -1) 42.59U (K -2) e(K)则,其中的参数有 耳=-0.7859, %=0.3679, b1=59.96, b2=42.59。本系统利用最小二乘递推算法(Recursive Least Squares, RL$进行辨识，将 辨识的结果与理论精确值进行比较，从而分析RLS算法的优劣。辨识模型选择N(K) =aN(K -1)-a2N(K -2) b1U (K -1) b2U (K -2) e(K)观测数据长度L = 480 ,则经最小二乘递推程序仿真后得到E0 al <1k478 double&g

10、t;45s47047147247347447547647747B4791|-0.7H7941-7879-07S79-0.7B79-0.7B7841-7876-D.7876-0.7E75-O.7S75-4).7874>田 a2 1x478 double*46947047147247347447547647747s479JjOlS&BSQ36SB03687Q.368&0.3&B5Q368603翊Q36B5-r田 bl <1x473 double>469470471472473474475476477478479-Al59.94:059.94&D5S

11、J45759.944259.933959.M5359.944259,9476我 05844田<bc47fl double>469470-4714724734744754764774784754工册W42.507442,5ff7742.5064_ 42,51W42.512142,513142.518743,517542.5J49图5估计的参数表1数据对比qa2b1实际数据-0.78590.367959.960042.5900估计数据-0.78490.388559.789942.7652706050aia2b1KOb2最小二乘递推算法辨识曲线40302010-10-20-30匚 050

12、100150200250300350400450500图6最小二乘递推算法辨识曲线附加代码如下：clc； clear;%对象参数a=1 -0.7859 0.3679;b=59.96 42.59;d=1;na=length(a)-1;nb=length(b)-1; %na、nb 为 A、B 阶次L=480; %仿真长度uk=zeros(d+nb,1);%输入初值：uk(i)表示 u(k-i)yk=zeros(na,1);%输出初值u=7*rand(L,1);%输入采用0-7v之间的随机数%u=7*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4,1);

13、 %输入幅值为 7V 的方波信号v=sqrt(100)*(rand(L,1)-0.5);%产生-0.50.5 之间的随机数 ，方差为 100theta=a(2) a(3) b(1) b(2);%对象参数真值thetae_1=zeros(na+nb+1,1); %thetaeu值P=10A6*eye(na+nb+1);for k=1:Lphi = -yk;uk(d:d+nb);%'止匕处 phi 为列向量y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+xi(k); % 采集输出数据 y(k)%递推最小二乘法K=P*phi/(1+phi'*P*phi);thetae

14、( : ,k)=thetae_1 + K*(y(k)-phi'*thetae_1);P=(eye(na+nb+1)-K*phi')*P;%更新数据thetae_1 = thetae( : ,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endplot(1:L,thetae); %line(1,L,theta,theta);xlabel('k');Iegend(,a1,a2',b1',b2,);title('最小二

15、乘递推算法辨识曲线,);2自适应控制部分2.1 提出控制问题直流电机在生活中、工业生产中等用的特别多，而电机的转动容易受输入电 压的影响，电压小小的波动就会带动电机转速的波动、超调等，针对这种情况， 文中在总结和分析前人的工作的基础上，提出将带有在线辨识的自适应算法应用 于电机控制系统中，旨在打破传统的控制策略，寻求更有效的控制方案将在线辨识自适应控制算法应用于该系统，利用系统的输入、输出信息，为了实现自适应控制和跟踪时变参数，采用递推最小二乘算法不断的循环调整PID参数值，使系统运行中保持合适的瞬态参数，以克服传统PID控制器参数不可改变的缺点，提高系统的控制品质，使直流电机系统具有更好的鲁

16、棒性和自适应能 力。2.2 自适应控制硬件直流电机控制系统硬件框图如图 2.1所示，其工作原理大致如下：Mcu作为 自适应控制的核心器件，得到经 A/D转换来的输入模拟信号，再结合编码器反 馈回来的直流电机的转速n(k),经MCU的自适应控制算法计算，输出信号 u(k) 给电子调速器，电子调速器可以将直流电压转换为三相电压信号送予直流电机， 其中，编码器是将电机的转速转换为脉冲信号反馈给Mcu o图2.1直流电机控制系统硬件框图2.3 算法推导及Matlab仿真系统采用自适应极点配置PID控制算法对直流电机的转速进行自适应控制， 控制算法结构框图如图2.2所示，控制系统由被控对象、控制器、辨识

17、器构成， 辨识器根据最小二乘递推算法，在线地计算被控对象的未知参数，控制器根据辨 识参数调整PID参数，经过不断的辨识和调整，使被控系统的性能指标渐近一 致地趋于最优。图2.2具有在线辨识的直流电机自适应控制系统被控对象为：N(K) =0.7859N(K -1) -0.3679N(K -2) 59.96U (K -1) 42.59U (K - 2) e(K) 式中e(K)为白噪声。下面首先首先实现极点配置PID算法，即离线控制，然后再这个基础之上， 实现在线系统辨识的自适应极点配置 PID控制算法。设被控对象为A(z)y(k) = z'B(z')u(k) +e(k)式中，u(k

18、)和y(k)表示系统的输入和输出，e(k)为外部扰动，d之1为纯延 时，且 112A(z ) = 1+a1za2z1_nB(z,)=b0 6z'dz， bnbz b对于本文的直流电机，d=1, nb =1。令期望闭环特征多项式为：A(z,):1式中，a和am2可以根据连续系统的特征多项式 &s) = s2+2,二6ns+0n2离散 化得到，即a2e y cos :111am2 = e为二其中，仃="蜀1 _，2,中=3工小- , ts为采样周期。根据期望系统模型参数计算期望闭环系统特征多项式，本系统中，令无阻尼 自然震荡角频率n = 5rad/s,系统阻尼比=0.70

19、7。由二 一 一；。1 - 2 =2.5一- -T -,1 -= 0.35am1 =-2e - cos=-0.0103a- e 2 - e -5a m2 一 e- e所以期望的闭环特征多项式为A(z,)=1 0.0103z，e2下面求自适应PID控制器：考虑如下反馈控制率呐二繇 yr(k) 一基 y(k)式中，yr (k)是被控过程参考输入信号。F(z)= f。 iz' "T且有,G(z)=1 uzgnkz”kH(z,)=h。几z也?箝当采用增量式PID控制器时，有如下公式成立，u(k) =Kpe(k) -e(k -1) K&k) Kde(k) -2e(k-1) e(

20、k - 2)式中，e(k) = yr (k) - y(k), >u(k) = u(k) -u(k -1)考虑上一时刻控制增量Au(k-1)的影响，使控制量更加平稳的增量式PID控制器，(1 : z).：u(k) =Kpy(k -1)-y(k)Kiy,(k)-y(k)Kd2y(k-1)-y(k) y(k-2)G(z') =1 (: -1)z， z'而且，F(z,)=(Kp K Kd) (-Kp -2Kd)zKdz萃H(z)=K在直流电机自适应PID控制系统中，期望电机转速为1000n/min ,输出电压 限幅15V。实现的算法如下：已知：模型阶次.,n及d,期望闭环特征多项

21、式；Step1设置初值P(0)、8(0)及遗忘因子九，输入初始数据和扰动；Step2采样当前输出y(k)和期望输出y.(k);Step3利用遗忘因子递推增广最小二乘法在线实时估计常，即A、B?;Step4求解单步Diophantine方程，得到多项式F, G, R的系数；Step5代入上式计算最小方差控制率u(k);Step6返回Step2,继续循环。根据以上步骤，首先给出对象参数已知时的极点配置 PID控制算法结果，如 图2.3所示。- - -O o o o o O 2 1-1-2 ,压电出输器制控J T,中俨 i卜弧军350O 50Q Q 5 3 O 5 2 2 150 00图2.3对象参

22、数已知时的仿真结果从图2.3可以看出，通过极点配置PID控制算法后的直流电机的转速可以很 好的跟踪给定信号，而控制器的输出也会迅速根据给定信号来调节相应的参数， 图中，因为加入了白噪声，所以有少许的波动。为了更好的调节直流电机，下面通过在线辨识系统参数，再使用自适应 PID 控制算法，其仿真结果如图2.4。图2.4对象参数未知时的仿真结果压昌出输器制控40002000速 转0-200001期望转速实际转速一 L 一-_rRU_50100150200250300350400450500步数观察图2.4,可以明显的观察到在前面几步控制中电机的输出转速n变动非常大，控制效果非常差，但是当这几步通过之

23、后，RLS算法辨识出直流电机的系 统参数后，经由PID控制算法，可以较好地控制电机的行为。比较图2.3和图2.4,在线辨识与离线辨识的效果明显不同，前者除了前几 步控制效果不太好之外，后期其优势明显强于后者，即在对象参数未知的情况下可以在线辨识系统并根据期望的闭环极点，配置 PID参数进行控制，系统适应 性非常强。附加代码：%PID极点配置控制(二阶系统，对象参数已知)clca=1 -0.7859 0.3679; b=59.96 42.59; d=1;Am=1 -0.0103 exp(-5); %对象参数及期望闭环特征多项式na=length(a)-1;nb=length(b)-1;nam=l

24、ength(Am)-1; %na、nb、nam 为多项式 A、B、Am 阶次 nf1=nb+d+2-(na+1)+1; ng=2; %nf1=nf+1%求解Diophantine方程，得到F、G、RF,G=diophantine(conv(a,1 -1),b,d,1,Am); %A0=1F1=conv(F,1 -1); R=sum(G);L=500; %控制步数uk=zeros(d+nb,1); %俞入初值：uk(i)表示 u(k-i)yk=zeros(na,1); %输出初值yr=1000*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4,1);

25、%期望输出，方波信号% e=2*ones(L,1); %常值干扰%加入白噪声告干扰%e2=2*ones(L,1);%labda=random('Normal',0,1,L,1)*20;%for i=1:L% e(i)=e2(i)*labda(i); %end% e=sqrt(1000)*randn(L,1);% 白噪声序列e=sqrt(0.1)*(rand(L,1)-0.5); %产生-0.50.5 之间的随机数 ，方差为 100 for k=1:Ltime(k)=k;y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+e(k); % 采集输出数据u(k)=(-F1

26、(2:nf1+1)*uk(1:nf1)+R*yr(k)-G*y(k);yk(1:ng)/F1(1); %求控制量%输出电压限幅if u(k)>=15u(k)=15;else if u(k)<=-15u(k)=-15;endend%更新数据for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endsubplot(2,1,1);plot(time,yr(1:L),'r:',time,y);xlabel('步数');ylabel('转

27、速');legend(期望转速,'实际转速');grid onsubplot(2,1,2);plot(time,u);xlabel('步数');ylabel('控制器输出(电压)')；grid on%自校正PID控制(二阶系统，对象参数未知)clca=1 -0.7859 0.3679; b=59.96 42.59; d=1;Am=1 -1.3205 0.4966; %对象参数及期望闭环特征多项式na=length(a)-1; nb=length(b)-1; nam=length(Am)-1; %na nb、nc、nam为多项式 A、B、C

28、、Am阶次nf1=nb+d+2-(na+1)+1; ng=2; %nf1=nf+1L=500; %控制步数uk=zeros(d+nb,1); %俞入初值：uk(i)表示 u(k-i)yk=zeros(na,1); %输出初值yr=1000*ones(L/4,1);-ones(L/4,1);ones(L/4,1);-ones(L/4,1); %期望输出%加入白噪声告干扰%e2=2*ones(L,1);%labda=random('Normal',0,1,L,1)*20;%for i=1:L% e(i)=e2(i)*labda(i);%end% e=sqrt(1000)*randn(L,1);% 白噪声序列e=sqrt(0.1)*(rand(L,1)-0.5);%产生-0.50.5 之间的随机数 ，方差为 100%RLS初值thetae_1=0.001*ones(na+nb+1,1); %初始化P=10A6*eye(na+nb+1);for k=1:Ltime(k)=k;y(k)=-a(2:na+1)*yk+b*uk(d:d+nb)+e(k); % 采集输出数据%递推最小二乘法phie=-yk(1:na);uk(d:d+nb);K=P*phie/(1+phie'*P*phie);thetae(:,k