一级直线倒立摆matlab程序

1、非线性作业一级直线倒立摆如图 1 所示 系统里的各参数变量M ：小车系统的等效质量（ 1.096kg）； m1 ：摆杆的质量（ 0.109kg）； m2 ：摆杆的半长（ 0.25m）；J：摆杆系统的转动惯量（ 0.0034kg*m）； g：重力加速度（ 9.8N/Kg ）； r：小车的水平位置（ m）；：摆角大小（以竖直向上为 0 起始位置，逆时针方向为正方向） ；Fh ：小车对摆杆水平方向作用力（ N）（向左为正方向）， Fh 是其反作用力； Fv ：小车对摆杆竖直方向作用力（ N）（向上为正方向）， Fv 是其反作用力； U：电动机经传动机构给小车的力，可理解为控制作用 u（向左为正方向）

2、； xp ：摆杆重心的水平位置（ m）； yp ：摆杆重心的竖直位置（ m）。11 一级倒立摆的数学建模 定义系统的状态为 r,r, , 经推导整理后可以达到倒立摆系统的牛顿力学模型：1）2）（I ml 2 ） mr l cos mgl sin2ml cos （M m）r mlsin u因为摆杆一般在工作在竖直向上的小领域内=0，可以在小范围近似处理：2cos 1,sin 0,0 ，则数学模型可以整理成：（I ml 2 ） mrl mgl（ 3）4）ml (M m)r u系统的状态空间模型为7)8)02g22mMm)(M(M gl m(M0010r0100r0r000.62930r 0.883

3、2+u000100027.828502.3566rr1000 r0u代人实际系统的参数后状态方程为：0010 01.2 滑模变结构在一级倒立摆系统的应用主要包括切换函数的设计、 控制率的设计和系统消除抖振的抑制。 基于线性二次 型最优化理论的切换函数设计，定义系统的优化积分指标是：J 0 x T Qxdt Q0，本文采用指数趋近律： S kS sgn(S) ，其中k和为正数。将其代人 S=Cx=0中，可以得到：S Cx CAx CBu kS sgn(S)( 9)控制率为： u (CB) 1(CAx kS sgn(S)(10)的选取主要是为了抑制系统的摩擦力和近似线性化所带来的误差和参数摄动等

4、因素，从而使得系统具有良好的鲁棒性。文中 k=25, =0。.8取变换矩阵 T。1000010(I ml2 )其中 T= 0ml00100001去 Q*11=diag(300 50 350)， Q22=10关于 Riccati方程 PA11 A11T P PA12Q221A1T2P Q11 0的解有 MATLAB 的 lqr函数可以解出692.7165 214.9281 359.6098P 214.9281 112.1894 149.8204359.6098 149.8204 354.4007C0100.0000 69.2717 210.6295 35.9610TS(x) C0 r r二 程序

5、% 主程序：直线一级倒立摆clear allclose allglobal C M0 Fts=0.02; %采样时间T=30; %仿真时间TimeSet=0:ts:T;para=;%options 为解微分函数中的调整参数 ,reltol 和 abstol 分别是设置相对误差和绝对 误差options=odeset(RelTol,1e-3,AbsTol,1e-3 1e-3 1e-3 1e-3); %options=;x0=0.5,0.3,0,0; %初始值 t,xout=ode45(daolibai2eq,TimeSet,x0,options,para);%固定格式，子程序调用 %返回值x1=

6、xout(:,1);x2=xout(:,2);x3=xout(:,3);x4=xout(:,4);s=C(1)*x1+C(2)*x2+C(3)*x3+C(4)*x4; % 切换函数或切换面 %选择控制器if F=1% for k=1:1:T/ts+1M0=40;%u(k)=-M0*sign(s(k);u=-M0*sign(s);%endelse if F=2beta=30;delta=0;for k=1:1:T/ts+1u(k)=-beta*(abs(x1(k)+abs(x2(k)+abs(x3(k)+abs(x4(k)+delta)*sign(s(k); endend%绘图figure(1)

7、;plot(t,x1,r);xlabel(time(s);ylabel(Cart Position);figure(2);plot(t,x2,r);xlabel(time(s);ylabel(Pendulum Angle);figure(3);plot(t,s,r);xlabel(time(s);ylabel(s);figure(4);plot(t,u,r);xlabel(time(s);ylabel(u);end% 子程序function dx=DxnamicModel(t,x,flag,para)% 自定义动态函数 global C M0 F%倒立摆经计算后模型 M=1.096;m=0.1

8、09;b=0.1;l=0.25;I=0.0034;T=0.005;g=0.98; k22=-(I+m*l2)*b/(I*(M+m)+M*m*l2);k23=m2*g*l2/(I*(M+m)+M*m*l2); k42=-m*l*b/(I*(M+m)+M*m*l2); k43=m*g*l*(M+m)/(I*(M+m)+M*m*l2);b12=(I+m*l2)/(I*(M+m)+M*m*l2); b14=m*l/(I*(M+m)+M*m*l2);A=0,1,0,0;0,k22,k23,0;0,0,0,1;0,k42,k43,0;b=0;b12;0;b14;%Ackermanns formula n1=-1;n2=-2;n3=-3;C=0,0,0,1*inv(b,A*b,A2*b,A3*b)*(A-n1*eye(4)*(A-n2*eye(4)*(A-n3*eye(4) );s=C*x;F=2;if F=1M0=40;u=-M0*sign(s);elseif F=2beta=30;delta=0;u=-beta*(abs(x(1)+abs(x(2)+abs(x(3)+abs(x(4)+delta)*sign(s);end%state equation dx=zeros(4,1);f0=0.5;ft=f0*sin(3*t);dx=A