倒立摆系统的建模及Matlab仿真

倒立摆系统的建模及Matlab 仿真1.系统的物理模型考虑如图(1)面内运动的二维问题。

图(1)倒立摆系统假定倒立摆系统的参数如下。

摆杆的质量：m=0.1g摆杆的长度：l =1m 小车的质量：M=1kg 重力加速度：g=9.8m/2s 摆杆的质量在摆杆的中心。

设计一个控制系统，使得当给定任意初始条件(由干扰引起)时，最大超调量δ ≤10%，调节时间ts ≤4s ，通过小车的水平运动使倒立摆保持在垂直位置。

2.系统的数学模型2.1建立倒置摆的运动方程并将其线性化。

为简化问题，在数学模型中首先假设：1)摆杆为刚体；2）忽略摆杆与支点之间的摩擦；3）忽略小车与接触面间的摩擦。

设小车瞬时位置为z,摆心瞬时位置为（θsin l z +）,在u 作用下，小车及摆均产生加速远动，根据牛顿第二定律，在水平直线远动方向的惯性力应与u 平衡，于是有u l z dtd m dt z d M =++)sin (2222θ 即： u ml ml z m M =-++θθθθsin cos )(2&&&&&① 绕摆轴转动的惯性力矩与重力矩平衡，因而有θθθsin cos )sin (22mgl l l z dt d m =⋅⎥⎦⎤⎢⎣⎡+ 即： θθθθθθθsin cos sin cos cos 22g l l z =-+&&&&&② 以上两个方程都是非线性方程，为求得解析解，需作线性化处理。

由于控制的目的是保持倒立摆直立，在试驾合适的外力条件下，假定θ很小，接近于零时合理的，则1cos ,sin ≈≈θθθ，且可忽略θθ2&项。

于是有u ml z m M =++θ&&&&)( ③ θθg l z =+&&&&④ 联立求解可得u Ml Ml m M u MM mg z 1)(1-+=+-=θθθ&&&&2.2列写系统的状态空间表达式。

选取系统变量4321,,,x x x x ， []Tx x x x x 4321,,,=则u Mlx Ml m M x x x u Mx M mg x x x 1)(134433221-+==+-==&&&&即[]Cxx x y Bu Ax u Ml M x Mlg m M Mmgz z dt d x ===+=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡-+⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡+-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=000110100)(0010000000011θθ&&& 代入数据计算得到：[][]0,0001,1010,01100100001000010==-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=D C B A T3.设计控制器3.1判断系统的能控性和稳定性[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----==01101110100101101032B A BA AB BQ k ，rank(k Q )=4,故被控对象完全可控 由特征方程 0)11(22=-=-λλλA I 解得特征值为 0，0，11±。

出现大于零的特征值，故被控对象不稳定3.2确定希望的极点希望的极点n=4，选其中一对为主导极点1s 和2s ，另一对为远极点，认为系统性能主要由主导极点决定，远极点只有微小影响。

根据二阶系统的关系式，先确定主导极点1.021≤=--ςπξσe p 可得59.0≥ξ，于是取6.0=ξ；取误差带02.0=∆有ns t ξω4=，则 1.67=n ω，闭环主导极点为22,11ξξω-±-=j s n =-1±0.8j,远极点选择使它和原点的距离大于主导极点与原点距离的5倍，取154,3-=s3.3采用状态反馈方法使系统稳定并配置极点状态反馈的控制规律为kx u -=，[]321k k k k k =；状态反馈系统的状态方程为Bv x BK A x+-=)(&，其特征多项式为 0122033141010)11()()(k k k k k k BK A I ----+-+=--λλλλλ ⑤希望特征多项式为3692.49964.28632)8.01)(8.01()15(2342++++=++-++λλλλλλλj j ⑥比较以上两式系数，解得状态反馈矩阵[]92.8154.33492.499.36----=k4.设计全维观测器4.1判断系统的能观性[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--==100001000010001)()(32C A C A C A CQ T T T g ，rank(g Q )=4,故被控对象完全可观 4.2确定观测器的反馈增益全维观测器的动态方程为GCx Bv x GC A x++-=)&))(；其特征多项式为 )11()11()11()(312021304g g g g g g GC A I --+--+-++=--λλλλλ ⑦取观测器的希望极点为：-45，-45，-3+3j ，-3-3j ；则希望特征多项式为3465013770258396)8.01)(8.01()15(2342++++=++-++λλλλλλλj j ⑧比较以上两式系数，解得观测器反馈矩阵[]TG 6498414826259496--=5.降维状态观测器的设计5.1建立倒置摆三维子系统动态方程设小车位移z 由输出传感器测量，因而无需估计，可以设计降维（三维）状态观测器，通过重新排列被控系统变量的次序，把需由降维状态观测器估计的状态变量与输出传感器测得的状态变量分离开。

将z 作为第四个状态变量，则被控系统的状态方程和输出方程变换为[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡z z y u z z z z dt d θθθθθθ&&&&&&1000010100010011001000010 ⑨简记为：[]⎥⎦⎤⎢⎣⎡==⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡211212122211211210x x I y y u b b x x A A A A x x &&式中[]T zx θθ&&=1，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=01110001011A ，[]T A 00012=，T b ]101[1-= y z x ==2，[]TA 00121=，21A =0，02=b ，11=I被控系统的n-q 维子系统动态方程的一般形式为v x A x +=1111&，121x A z ='式中u b u b y A v 1121=+=,z y u b y A y z &&&==--='222z '为子系统输出量。

故倒置摆三维子系统动态方程为[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡='⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-+⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡θθθθθθ&&&&&&z z u z z dt d 00110101101000105.2.判断子系统的可观测性A1=[0 -1 0;0 0 1;0 11 0];C1= [1 0 0];Qg1=obsv(A1,C1);r=rank(Qg1)运行Matlab 程序；结果为r=3,故该子系统可观测 降维状态观测器动态方程的一般形式为()()()[]yh x y A h A h A h A u b h b A h A +=-+-+-+-=ωωω122122*********&&式中h=[]Th h h 210。

考虑被控对象参数，单倒置摆降维观测器动态方程的一般形式为y h h h x y h h h h h h h h u h h h ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡+=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡+-+---+⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-+⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----=2101120210120210111010111001ωωω&&^`5.3确定三维状态观测器的反馈矩阵h三维状态观测器的特征多项式为()()()20120321111111h h h h A h A I --+--++=--λλλλ设希望的观测器闭环极点为-45，-3+3j ，-3-3j ，则希望特征多项式为()()()8102885133334523+++=-++++λλλλλλj j比较以上两式系数，解得h=[]1371-29951-故所求三维状态观测器的动态方程为y y x x y u ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-+⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-+⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=113712995100000100001000016663213878230210101113711029901511ωωω&)& 6.Matlab 仿真分析6.1源程序通过Matlab 对用全维状态观测器实现状态反馈的倒置摆系统进行仿真分析，下面是文件名为Inversion_pendulum_system.m 的源程序 %倒立摆系统建模分析%a)判断系统能控性和能观性 clear all ; clcA=[0 1 0 0;0 0 -1 0;0 0 0 1;0 0 11 0]; B=[0;1;0;-1]; C=[1 0 0 0]; D=0;Uc=ctrb(A,B);rc=rank(Uc); n=size(A); if rc==ndisp('The system is controlled.') elseif rc<ndisp('The system is uncontrolled.') endVo=obsv(A,C); ro=rank(Vo); if ro==ndisp('The system is observable.')^` elseif ro~=ndisp('The system is no observable.')end%b)判断系统稳定性P=poly(A),v=roots(P)Re=real(v);if(length(find(Re>0))~=0)disp('The system is unstable and the ubstable poles are:')v(find(Re>0))elsedisp('The system is stable!');end% c)极点配置与控制器-全维状态观测器设计与仿真pc=[-1+0.8*j,-1-0.8*j,-15,-15];po=[-45 -45 -3+3*j -3-3*j];K=acker(A,B,pc),G=acker(A',C',po)'Gp=ss(A,B,C,D); %将受控过程创建为一个LTI对象disp('受控对象的传递函数模型：');H=tf(Gp)Af=A-B*K-G*C;disp('观测器——控制器模型：');Gc=ss(Af,-G,-K,0) %将观测器-控制器创建为一个LTI对象disp('观测器——控制器的极点：');f_poles=pole(Gc)GpGc=Gp*Gc; %控制器和对象串联disp('观测器——控制器与对象串联构成的闭环系统模型：');Gcl=feedback(GpGc,1,-1) %闭环系统disp('闭环系统的极点和零点：');c_poles=pole(Gcl)c_zeros=tzero(Gcl)lfg=dcgain(Gcl) %低频增益N=1/lfg % 归一化常数T=N*Gcl; %将N与闭环系统传递函数串联x0=[100 10 30 10 0 0 0 0];%初始条件向量t=[0:0.01:1]'; %时间列向量r=0*t; %零参考输入[y t x]=lsim(T,r,t,x0); %初始条件仿真plot(t,x(:,1:4),'-.',t,x(:,5:8)) %由初始条件引起的状态响应title('\bf状态响应')legend('x1','x2','x3','x4','x1hat','x2hat','x3hat','x4hat')figure(2)step(T)title('\bf阶跃响应')figure(3)impulse(T)title('\bf脉冲响应')^` 6.2 程序运行结果The system is controlled.The system is observable.P =1 0 -11 0 0v =3.3166-3.3166The system is unstable and the ubstable poles are:ans =3.3166K =-36.9000 -49.9200 -334.5400 -81.9200G =962594-14826-64984受控对象的传递函数模型Transfer function:s^2 - 1.776e-015 s - 10-----------------------s^4 - 11 s^2观测器——控制器模型：a =x1 x2 x3 x4x1 -96 1 0 0x2 -2557 49.92 333.5 81.92x3 1.483e+004 0 0 1x4 6.495e+004 -49.92 -323.5 -81.92b =u1x1 -96x2 -2594x3 1.483e+004x4 6.498e+004c =^` x1 x2 x3 x4y1 36.9 49.92 334.5 81.92d =u1y1 0Continuous-time model.观测器——控制器的极点：f_poles =1.0e+002 *-1.4948 + 1.8786i-1.4948 - 1.8786i1.7424-0.0328观测器——控制器与对象串联构成的闭环系统模型：a =x1 x2 x3 x4 x5 x1 0 1 0 0 0x2 0 0 -1 0 36.9x3 0 0 0 1 0x4 0 0 11 0 -36.9x5 96 0 0 0 -96x6 2594 0 0 0 -2557x7 -1.483e+004 0 0 0 1.483e+004x8 -6.498e+004 0 0 0 6.495e+004x6 x7 x8x1 0 0 0x2 49.92 334.5 81.92x3 0 0 0x4 -49.92 -334.5 -81.92x5 1 0 0x6 49.92 333.5 81.92x7 0 0 1x8 -49.92 -323.5 -81.92b =u1x1 0x2 0x3 0x4 0x5 -96x6 -2594x7 1.483e+004x8 6.498e+004c =^`x1 x2 x3 x4 x5 x6 x7 x8 y1 1 0 0 0 0 0 0 0 d =u1 y1 0Continuous-time model. 闭环系统的极点和零点： c_poles =-45.0000 -45.0000 -15.0001 -14.9999 -3.0000 + 3.0000i -3.0000 - 3.0000i -1.0000 + 0.8000i -1.0000 - 0.8000i c_zeros = 3.1623 0.2263 -0.1707 -3.4312 -3.1623 lfg =1.0000 N =1.0000由控制器——全维状态观测器实现的倒立摆系统在初始条件下引起的状态变量的响应、输出变量的阶跃响应和脉冲响应如下图图(2)状态响应()t x （虚线）和()t x )（实线）^` 图(3)阶跃响应()t y图(4)脉冲响应()t y。