系统建模

题目： 假设：

1. 自行车车轮的转动惯量，质量和宽度忽略不计。滚动过程无侧向和横向滑动； 2. 自行车平衡时以匀速率前进，车体在前进中遇到小扰动而偏离平衡位置； 3. 将连接前后轮的三脚架视为其质量集中在其质心，同时把前后轮胎视为质心

在车轮圆心的圆环；

4. 自行车在水平面内运动，后轮的滚动速度为常量。

行驶中的自行车可以近似为倒立摆模型，自行车在运行过程中的各视图如图1和图2所示。

Fmg 图1 车体正视图

图2 车体右侧视图和俯视图

图中各符号参数及相应值如下：

表1 各符号参数及相应值

物理意义 车体质量 后轮着地点到车体质心距离 前轮着地点到车体质心距离 质心高度 车把转动惯量 车体转动惯量 受扰动后车体偏离垂直方向角度 后轮速度 转弯速度 自行车转弯轨迹中心 表示符号 m 值 18kg 0.81m 1.17m 0.7m 0.153kg•m2 12.0 kg•m2 b a h Ih I  v0 v o 车后轮转弯时的曲率半径 质心运动路径瞬时曲率半径 车体质心的夹角 前轮中心的夹角 输入量 输出量 r0 r     假定无侧滑由俯视图车体的运行轨迹由前后距离和车把的转角决定

br0tanrsinar0tan

  是可控制量

v0r0v0atanvrddt

实际情况中前后轮及质心速度不同，此处假定三者相等，由于瞬间曲率半径远大于前后轮着地点距离a，r0r，质心速度v满足

vrddtbv0tanasinba

vyvsinv0tan车体转弯时，作用于质心的离心力

Fmv2rcosm2dvydtbv02m(v0tanad

)acosdt重力分量：mgsin

ddt22ImghsinFhcos

问题：

1) 建立自行车运动方程，输入量为，输出量为。

2) 在建好的模型基础上，建立状态空间方程，x1，x2。利用MATLAB

画出向量场，指出系统的平衡点。

3) 在平衡点附近，近似得到线性化之后的模型，写出系统的传递函数。 4) 利用李雅普诺夫第二方法分析系统的平衡点特性。

5) 利用离散化公式，将系统离散化并化为最小二乘形式rH。 6) 假定路面具有未知摩擦力f，设计此时的最小二乘格式。 7)

f5时，用

MATLAB编程，辨识得到f的估计值，递推算法实现。

ck(BkA)0，分析稳定性条8) 假定搬动车把的控制率为k，件，求出维持系统平衡的最小速度v0，解释为什么较快速度下易维持平衡。

（1）建立模型

根据假设条件，对自行车车体进行受力分析，可以得到自行车受到的离心力F,

Fm(v0tana2bv02dacosdt) (1.1)

取与偏离方向相反的方向为正方向，得运动方程如下

Iddt22Fhcosmghsin (1.2)

将(1.1) 带入 (1.2) 有

v0tana2bv02dacosdtgtanId22mhcosdt (1.3)

将b=0.81,a=1.17,Ih=0.153,I=12.0,m=18kg,h=0.7带入（1.3）得

0.85v0tan0.692v02dcosdt10tan0.96dcosdt22 (1.4)

上式即为自行车的运动方程，其中输入量为，输出量为。

（2）确定平衡点

在数学模型基础上，建立状态空间方程如下

.x1x2 (2.1) .v0d210.42sinx1x20.cosx1v0tan0.72cosx12cosdt

当输入量为零时，状态方程可化简如下

.xx21..x210.42sinx1 (2.2)

利用matlab仿真向量场 程序代码为：

x0=-3.6:.5:3.7; y0=-3.7:.5:3.7;

[x1,x2]=meshgrid(x0,y0);

d=sqrt((x2.^2)+(-10.42*sin(x1)).^2); u=x2./d;

v=(-10.42*sin(x1))./d; hold on

quiver(x1,x2,u,v,0.5,'b');

向量场仿真结果为:

图1

系统的平衡点为漩涡的中心，共有3处。根据方程求解可知系统的平衡点分别为（0,0）、0，由实际情况可知，自行车平衡点只能取0，即平衡点只有一个(0,0)。

（3）在平衡点处，因为

由

cos1,sin,tan,cos1 得到

ddt22mghI总mhv20I总amhbvI总a0ddt （3.1）

进行拉普拉斯变换得传递函数为：

G(s)mhbvsmhv002 （3.2）

(s)(s)2I总asmghI总（4）构造正定函数

22 V1(x)x1x2 (4.1)

V1(x)2x ) (4.2) x22x2x22x(1x10.42 s 1xi n11在平衡点处，近似等于18.84x1x2

...

，故当x1，x2异号时V(x)为正，系统不稳定;

.1当他们同号时.V1(x)为负,系统稳定。

（5）因为

dd22mv0h0dtdt2mghmhvIaIaI 所以离散化后得

(k+2)2(k+1)(k)2(k)]T211.06(k)0.v0(k)0.72v0[(k+1)T 取T=0.2可进一步得

(k+2)2(k+1)0.5576(k)+(0.0356v20-0.144v0)(k)0.144v0(k+1) (k)2(k1)0.5576(k2)+(0.0356v20-0.144v0)(k2)0.144v0(k-1)令

H[(k1),k(2),k(1k), ( [2,0.5576,v020.144v35vT0,0.00 0 . 1 4 4 ] 则

(k)H （6）当有常摩擦力f时，其运动方程如下

2Iddt2Fhcosmghsinfhcos 将F带入（6.1）有

22mv0tandambv0acos2dtmgtanIdhcosdt2f 将其线性化，得

2mv0bv0dId2amadtmghdt2f 带入相关数据，则(6.3)可化简为

（5.1）

5.2）

5.3）

（5.4）

5.5）

5.6）

5.7）

(6.1)

(6.2)

(6.3)

（ （ （（（

0.085v0.069v020ddt0.096ddt22f (6.4)

再将其离散化，可得

(k)0.096(k2)0.085v0(k)0.069v0(k1)0.006f (6.5)

2进而，有

(nk)0.096(nk2)0.085v0(nk)0.069v0(nk1)0.006fe(nk)2(6.6)

令

k[(nk2),(nk),(nk1),1] (6.7) TT[0.096,0.085v0,0.069v0,0.006f]

2 (6.8)

则

1T2TNTN(n1)(n)(nN2)(n1)(n2)(nN1)1(n1)1 (6.9) (nN1)1(n)yNy(n1)y(n2) (6.10)

y(nN)则此时的最小二乘格式为

T1TN(NN)NyN (6.11) 

（7）用递归最小二乘法辨识得到f估计值，程序流程图如下

θ(k)给被辨识参数θ和p赋初值输入n组初始数据计算K（k）分离参数计算θ(k)

画出被辨识参数的各次递推估计值的图形计算P（k）画出被辨识参数的相对误差的图计算被辨识参数的相对变化量停N参数收敛满足要求Y图2

令f5用matlab编程辨识得到f的估计值（递推方法实现） 已知f5N 代入数据得：

(k)(k1),(k2),(k1),(k2),0.00922,0.6257,1.4355,2.109,5

T递推算法的matlab的程序见附录 运行结果如下：

0.03 0.02 0.01 0 -0.01 -0.02 -0.030

Parameter Identification with Recursive Least Squares Method543Identification Precision5000-500210-1-2000-2-3-2500051015-1000-150051015051015 输入序列图 递推最小二乘参数辨识图 辨识精度图

（8）由题知

vrkddt (8.1)

所以

vrk进而可得

vCkv(BkA)v0

....ddt22 (8.2)

(8.3)

进行拉氏变换得

sv(s)Cksv(s)(BkA)v(s)02 (8.4)

由劳斯判据可知

kAB (8.5)

ddtAdBdt所以

ddtkddt,vrk,v (8.6)

即维持系统平衡的最小速度为

AdBdt

又由F知，当v越大时，离心力越大，因此以较快速度行驶易维持平衡。

附录1

递推算法辨识f的程序

clear %清理工作间变量 L=15; % M序列的周期

f=5;

y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值 for i=1:L; %开始循环，长度为L

x1=xor(y3,y4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”

x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出 x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出 x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出

y(i)=y4; %取出第四个移位寄存器的幅值为\"0\"和\"1\"的输出信号，即M序列 if y(i)>0.5,u(i)=-0.03; %如果M序列的值为\"1\辨识的输入信号取“-0.03” else u(i)=0.03; %如果M序列的值为\"0\辨识的输入信号取“0.03” end %小循环结束

y1=x1;y2=x2;y3=x3;y4=x4; %为下一次的输入信号做准备 end %大循环结束，产生输入信号u figure(1); %第一个图形

stem(u),grid on %显示出输入信号径线图并给图形加上网格 z(2)=0;z(1)=0; %设z的前两个初始值为零 %v=randn(15,1); for k=3:15; %循环变量从3到15

z(k)=2*z(k-1)+0.6257*z(k-2)+1.4355*u(k-1)+2.109*u(k-2)-0.0092*f; %输出采样信号 end

%RLS递推最小二乘辨识

c0=[0.001 0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量

p0=10^6*eye(5,5); %直接给出初始状态P0，即一个充分大的实数单位矩阵 E=0.000000005; %取相对误差E=0.000000005 c=[c0,zeros(5,14)]; %被辨识参数矩阵的初始值及大小 e=zeros(5,15); %相对误差的初始值及大小 for k=3:15; %开始求K h1=[-z(k-1),-z(k-2),u(k-1),u(k-2) -0.0092]';

x=h1'*p0*h1+1;

x1=inv(x); %开始求K(k)

k1=p0*h1*x1; %求出K的值 d1=z(k)-h1'*c0; c1=c0+k1*d1; %求被辨识参数c

e1=c1-c0; %求参数当前值与上一次的值的差值

e2=e1./c0; %求参数的相对变化

e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1; %新获得的参数作为下一次递推的旧参数

c(:,k)=c1; %把辨识参数c 列向量加入辨识参数矩阵的最后一列 p1=p0-k1*k1'*[h1'*p0*h1+1]; %求出 p(k)的值 p0=p1; %给下次用

if e2<=E break; %如果参数收敛情况满足要求，终止计算 end %小循环结束 end %大循环结束

c,e %显示被辨识参数及其误差(收敛)情况 %分离参数 a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); ff=c(5,:) ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:);

eff=e(5,:);

figure(2); %第二个图形

i=1:15; %横坐标从1到15

plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':' ,i,ff,'b') %画出a1，a2，b1，b2的各次辨识结果 title('Parameter Identification with Recursive Least Squares Method') %图形标题 figure(3); %第三个图形

i=1:15; %横坐标从1到15

plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:',i,eff,'b') %画出a1，a2，b1，b2的各次辨识结果的收敛情况

title('Identification Precision') %图形标题