24小时热门版块排行榜

返回列表

lslhcqy

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 3.5
帖子: 8
在线: 1.6小时
虫号: 3720573
注册: 2015-03-09
性别: GG

[求助] 新手求教如何用bvp4c求解free flying robot这个问题

新手上路，求老司机带带啊！！！
问题描述如附图所示，我运行程序总是报错：HW_5bvpError using bvp4c (line 251)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in HW_5bvp (line 4)
sol=bvp4c(@HW_5ode,@HW_5bc,solinit,options);

另外，问题里边有个条件为Ti≤1，我不知道该咋用，程序里没有体现出来。搞不清楚它是状态变量，还是控制变量:'(

下边是我编写的程序，包括5个.m文件
1.function  SetGlobalu( v )                               %该函数是将控制变量u（u1,u2,u3,u4）在第5个函数中设置为全局变量
global u
u=v;
end

2.function r = GetGlobalu                               %该函数是为了在其它函数中共享第5个函数设置的控制变量u
global u
r=u;
end
                                                                     %该函数是程序的主函数
3.HW_5bvp.m
function HW_5bvp
solinit=bvpinit(linspace(0,1,10),[-10;-10;0;0;pi/2;0;0;1;0;2;0;1;-1;1],12); %设定时间初值及6个状态变量及其对应Px的初值，未知参数的初值设为12
options=bvpset('Stats','on','RelTol',1e-12,'AbsTol',1e-12);
sol=bvp4c(@HW_5ode,@HW_5bc,solinit,options);
y=sol.y;
time=sol.parameters*sol.x;
plot(time,y);
end

4.function res = HW_5bc(ya,yb,T )                      %设定问题的边界，y(1~13)分别表示x,y,Vx,Vy, θ, ω,Px,Py,PVx,PVy,Pθ,Pω，H（tf）
%哈密顿算子为H= Px* Vx+ Py*Vy+ PVx*cosθ*(u1-u2+u3-u4)+ PVy*sinθ*(u1-u2+u3-u4)+ Pθ*ω+ Pω*[α*(u1-u2)+ β*(u3-u4)]+(u1+u2+u3+u4)
r=GetGlobalu;
T3=r(1)-r(2)+r(3)-r(4);
T4=r(1)+r(2)+r(3)+r(4);
res=[ya(1)+10
ya(2)+10
ya(3)
ya(4)
ya(5)-pi/2
ya(6)
yb(1)
yb(2)
yb(3)
yb(4)
yb(5)
yb(6)
yb(7)*yb(3)+yb(8)*yb(4)+yb(11)*yb(6)+yb(12)*0.2*T3+yb(9)*cos(yb(5))*T3+yb(10)*sin(yb(5))*T3+T4];
                                                                                                            %y(13)=0,该值为横截条件，H（tf）=0，因为tf自由
end

5.function dydt = HW_5ode( x,y,T )
r=zeros(4,1);
if((1+y(9)*cos(y(5))+y(10)*sin(y(5))+y(12)*0.2)>0)       %这个是根据附图2里边的那个λTg来判断出来u=0或者1000的。λ即为前边H中的Px。
r(1)=0;r(3)=0;
end
if((1+y(10)*sin(y(5))+y(9)*cos(y(5))+y(12)*0.2)<0)
r(1)=1000;r(3)=1000;
end

if((1-y(9)*cos(y(5))-y(10)*sin(y(5))-y(12)*0.2)>0)
r(2)=0;r(4)=0;
end
if((1-y(9)*cos(y(5))-y(10)*sin(y(5))-y(12)*0.2)<0)
r(2)=1000;r(4)=1000;
end

SetGlobalu(r);
T1=r(1)-r(2);
T2=r(3)-r(4);
T3=r(1)-r(2)+r(3)-r(4);
dydt=T*[y(3)
y(4)
(T1+T2)*cos(y(5))
(T1+T2)*sin(y(5))
y(6)
0.2*(T1+T2)
0
0
y(7)
y(8)
-y(9)*sin(y(5))*T3+y(10)*cos(y(5))*T3
0];
end

问题描述.png

公式推导.png

回复此楼

» 猜你喜欢

你我生活在不同的世界

1楼 2016-12-01 08:38:33

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

相关版块跳转我要订阅楼主 lslhcqy 的主题更新

返回列表