本帖最后由 black1man 于 2018-11-29 23:57 编辑

我先说下我的仿真：

是对一个船舶系统进行反步滑模控制，其中simulink仿真中的闭环控制可以完美收敛，可是同样用m文件来做，却是彻底的发散。

(需要说一下，开环的情况下，simulink和m文件的仿真结果完全一致，只有在闭环情况下，一个收敛，一个发散)

请问：是我的m文件哪里编的不对吗？为什么会出现这种情况。程序和simulink都传到了附件里。

研究了一个晚上，完全没有发现问题所在，脑子都懵了，不知道哪里犯了错误，希望有朋友指点迷津。谢谢！

我把m文件贴了出来贴出来：(附件里有simulink和m程序，可以下载后直接运行)

clear;clc;format long;

T = 0.5; count = 200; step_T = 2;    %周期，周期数，周期步数

h = T/step_T;

%% 船体数据

M = 1e9* [ 0.0068 0 0; 0 0.0113 -0.0340 ; 0 -0.0340 4.4524];

D =  1e8* [0.0008 0 0 ; 0 0.0025 -0.0203 ; 0 -0.0203 3.8481];

L1=diag([0.22;0.22;22]); L2=diag([4e6;10e6;900e6]);   % 反步法的增益

omega_bar=[5e6;5e6;5e6];                                              % 干扰的上限

%% 变量初始化

eta=[10;10;1];

eta = J_func(eta);      %弧度转换到[-pi,pi]之间

nu= [0;0;0];

omega=[0;0;0];         %干扰

tau=[0;0;0];               %控制力

etad=[0;0;0];             %期望值

detad=[0;0;0];           %期望值的一次导

ddetad=[0;0;0];         %期望值的二次导

%% 画图

t_draw = zeros(1,count);

eta_draw = zeros(3,count);

omega_draw=zeros(3,count);

%% 主循环

for i = 1:1:count                      %主循环

t = (i-1)*T;

t_draw(i)    = t;              %       eta = J_func(eta);      %弧度转换到[-pi,pi]之间

eta_draw(:,i) = eta;

%% 干扰

omega = 4.85e6*[1+0.01*sin(0.01*t)+0.02*cos(0.2*t)*sin(0.01*t);

1+0.01*sin(0.01*t)+0.02*cos(0.2*t)*sin(0.01*t);

1+0.01*sin(0.01*t)+0.02*cos(0.2*t)*sin(0.01*t)];

omega_draw(:,i)=omega;

%% 控制

tau= tau_func(L1,L2,eta,nu,etad,detad,ddetad,M,D,omega_bar);

%% 船模

for j = 1:1:step_T                 %次循环 四阶龙格库塔

keta1= deta(eta,nu);

knu1 = dnu(M,D,nu,omega,tau);

keta2= deta(eta+h/2*keta1,nu+h/2*knu1);

knu2 = dnu(M,D,nu+h/2*knu1,omega,tau);

keta3= deta(eta+h/2*keta2,nu+h/2*knu2);

knu3 = dnu(M,D,nu+h/2*knu2,omega,tau);

keta4= deta(eta+h*keta3,nu+h*knu3);

knu4 = dnu(M,D,nu+h*knu3,omega,tau);

eta = eta + h/6*(keta1 + 2*keta2 + 2*keta3 + keta4);

nu  = nu +  h/6*(knu1 + 2*knu2 + 2*knu3 + knu4);

end

end

%% 画图

figure(1)

subplot(3,1,1)

plot(t_draw,eta_draw(1,:),'LineWidth',2);

subplot(3,1,2)

plot(t_draw,eta_draw(2,:),'LineWidth',2);

subplot(3,1,3)

plot(t_draw,eta_draw(3,:),'LineWidth',2);

figure(3)

plot(eta_draw(1,:),eta_draw(2,:),'LineWidth',2);

figure(2)             % 干扰

subplot(3,1,1)

plot(t_draw,omega_draw(1,:),'LineWidth',2);

subplot(3,1,2)

plot(t_draw,omega_draw(2,:),'LineWidth',2);

subplot(3,1,3)

plot(t_draw,omega_draw(3,:),'LineWidth',2);

%——————————————————————————

%以下是子函数

%——————————————————————————

function A = J_func(eta)%将[pi,3pi]和[-3pi,-pi]之间的转到-pi到pi之间

if eta(3)>pi

eta(3) = eta(3)-2*pi;

else if eta(3)

eta(3) = eta(3)+2*pi;

end

end

A = eta;

%————————————————————

function A =dnu(M,D,nu,omega,tau)

A = M\tau - M\D*nu + M\omega;

%————————————————————

function A = deta(eta,nu)

psi=eta(3);

J=[cos(psi),-sin(psi),0;

sin(psi),cos(psi),0;

0,0,1];

A=J*nu;

%————————————————————

function A =  tau_func(L1,L2,eta,nu,etad,detad,ddetad,M,D,omega_bar)

psi=eta(3);

J=[cos(psi),-sin(psi),0;

sin(psi),cos(psi),0;

0,0,1];

dJT= [-sin(psi), cos(psi),0;

-cos(psi),-sin(psi),0;

0,0,0];

z1=eta-etad;

dz1=J*nu-detad;

nud=J'*(detad-L1*z1);

dnud= dJT*(detad-L1*z1)+J'*(ddetad-L1*dz1);

z2=nu-nud;

s=z2;

A= D*nu +M*dnud-diag(sign(s))*omega_bar-L2*s-J'*z1;

untitled.jpg

(46.38 KB, 下载次数: 0)

2018-11-29 23:30 上传

simulink仿真闭环系统

untitled1.jpg

(50.65 KB, 下载次数: 0)

2018-11-29 23:33 上传

m文件闭环仿真

untitled2.jpg

(57.09 KB, 下载次数: 0)

2018-11-29 23:34 上传

m文件开环仿真

untitled3.jpg

(64.61 KB, 下载次数: 0)

2018-11-29 23:36 上传

simulink的开环仿真

2018-11-29 23:27 上传

点击文件名下载附件

2.21 KB, 下载次数: 17

2018-11-29 23:27 上传

点击文件名下载附件

1.38 MB, 下载次数: 26