ODE45 error that the returning vector has not the same length with the initial conditions.

Hi

Is there anyone can help me? When I call for the function file. It always came up like this message. I really do not know how to figure it out. Is there anyone can debug it? Thanks.

HW_4_MODEL returns a vector of length 3, but the length of initial conditions vector is 6. The vector returned by HW_4_MODEL and the initial
conditions vector must have the same number of elements.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> hw_4_animation at 71
[t,x] = ode45('hw_4_model',[0 6],x0);

function xdot = hw_4_model(t,x)
%
% Checked 3/8/11 rgl
%
global m Iz Cf Cr steer_angle L1 L2 Wf Wr mu A P_max Fd

Vx=x(1);Vy=x(2);omegaz=x(3);
Xt = x(4); Yt = x(5); Phit = x(6); 

if (t>=10 && t<=12)
    Fyfd = 0;
else
    Fyfd = 0;
end

% double lane change
if (t<1) deltaf = 0; end;
if (t>=1 & t<2) deltaf = steer_angle; end;
if (t>=2 & t<3) deltaf = -steer_angle; end;
if (t>=3 & t<4) deltaf = 0; end;
if (t>=4 & t<5) deltaf = -steer_angle; end;
if (t>=5 & t<6) deltaf = steer_angle; end;
if (t>=6) deltaf = 0; end;

alphaf = deltaf - atan((L1*omegaz+Vy)./Vx);
alphar = atan((L2*omegaz-Vy)./Vx);

Fyf = 2*Cf*alphaf;
Fyr = 2*Cr*alphar;

Fxf_max=mu*Wf;               % front traction limit
Fxr_max=mu*Wr;               % rear traction limit

P_max_f=P_max*Wf/(Wf+Wr);    % max front power
P_max_r=P_max*Wr/(Wf+Wr);    % max rear power

Fxf=min(Fxf_max,P_max_f./Vx*1000);     % front drive limited by both of traction and power
Fxr=min(Fxr_max,P_max_r./Vx*1000);     % rear drive limited by both of traction and power
Fd=Fxf+Fxr;                            % total drive force

Vx_dot=(m*Vy.*omegaz+Fxf.*cos(deltaf)+Fxr-Fyf.*sin(deltaf))/m;

Vy_dot = (-m*Vx.*omegaz + Fyr + Fyf.*cos(deltaf) + Fxf.*sin(deltaf) + Fyfd)/m;

omegaz_dot = (L1*Fyf.*cos(deltaf) - L2*Fyr + L1*Fxf.*sin(deltaf) + Fyfd*L1)/Iz;

Xtdot = Vx.*cos(Phit)-Vy.*sin(Phit);

Ytdot = Vx.*sin(Phit)+Vy.*cos(Phit);

phidot = omegaz;

xdot=[Vx_dot;Vy_dot;omegaz_dot;Xtdot;Ytdot;phidot];

%%%%%%%above is the function file, below is the M file call for %%%%%%%the function file.

clc
clear all
global m Iz Cf Cr B R_w steer_angle L1 L2 Fyfd tdon tdoff Wf Wr A P_max T_max Fd
% Baseline values
g = 9.81;
L = 3.075; L1 = 1.568; L2 = L-L1;
B=1.164;
R_w=13/39.37;
mu=0.8;
A=2.32;             % cross-section area, m^2
P_max=110;          % max power, kW
T_max=1000;         % max torque, N-m

% Inertia parameters
m = 1945;           % total mass, kg
W = m*g;      % weight, N

iyaw = 0.992;    % yaw dynamic index
% the following is a defined relation between yaw dynamic index and 
% the yaw moment of inertia (see Dixon reference and table of car specs)
Iz = iyaw*m*L1*L2;    % moment of inertia about z, kg-m^2

Wf = L2*W/L;    % static weight on front axle
Wr = L1*W/L;    % static weight on rear axle
% Refer to Wong, Section 1.4 for guide to the following parameters
CCf = 0.171*180/pi;  % front corning stiffness coefficient, /rad
CCr = 0.181*180/pi;  % rear corning stiffness coefficient, /rad

Cf = CCf*Wf/2;    % corning stiffness per tire, N/rad (front)
Cr = CCr*Wr/2;    % rear cornering stiffness per tire, N/rad (rear)

% Criteria for stability (Rocard basic model - no steer angle)
Rfactor = Cr*L2 - Cf*L1;
% critical velocity
% if Rfactor is positive, there is no critical speed
if Rfactor<0,
    Vc = sqrt(2*Cf*Cr*L*L/(m*(Cf*L1-Cr*L2)));  % critical speed
else
    Vc = 0;
end

% Steering wheel angle
steer_wheel = 50;       % degrees
steer_gain = 20;       % this was introduced to compare with examples in CarSim
steer_angle = steer_wheel*pi/180/steer_gain;

% Expected turn radius
if steer_angle==0
    sprintf('No steer angle')
    Fyfd = 1*100; % disturbance force in lateral direction on front tire
else
    R = L/steer_angle;
    Fyfd = 0;
end

% Initial conditions for states
Vxo=1*10^(-8);
Vyo = 0.0;
omegazo = 0;
Xt0 = 0;
Yt0 = 0;
Phit0 = 0;

x0=[Vxo;Vyo;omegazo;Xt0;Yt0;Phit0];

tend = 15;
tdon = tend/3;
tdoff = tdon*1.15;

[t,x] = ode45('hw_4_model',[0 6],x0);

figure(2)
subplot(1,2,2),plot(t,x(:,1)), title ('Longditudinal Velocity')

Answer 1

[Walter Roberson]

I would put in a breakpoint at the xdot= line, execute it, and see what size(xdot) came out to be. If it was not 6 x 1, I would start looking through the variables that make up xdot to see if any of them somehow ended up empty, and if so, I would backtrack from there to figure out how they ended up empty.

[Mingrui Wang]

actually, I am sure that something is wrong in the ODE below.

Vx_dot=(m*Vy.*omegaz+Fxf.*cos(deltaf)+Fxr-Fyf.*sin(deltaf))/m;

because this ODE was added for the new homework.The previous homework is about the rest 5 ODEs, and it works well. This time I need to add a new ODE. Also I have did some modifications to match the new 6-ODE file. However, it always come out the error message that I cannot move on. It ruined my spring-break.

Answer 2

[Jan]

A remark:

The ODE45 integrator is designed to integrate a smooth function. Your function has discontinuities caused by the non-smooth deltaf. This is a brute misusage of the intergrator and ODE45, which evaluates the function at different intermediate steps to calculate the average derivative for each full step. For the stepsize controller this means, that e.g. half of the intermediate values can belong to the t<1 interval and the other half to the t>=1 interval. The effects are not predictable and therefore it is a scientific sin to use ODE45 for discontinuos functions.

Surprisingly the results are not too bad, but why? ODE45 reduces the stepsize until the local discretization error "hides" the discontinuity due to cancellation errors. If the system is stable, the result can be near to the expected solution - accidently.

A scientific and serious simulation uses the event function to stop the integration at any discontinuity and restart the integrator. Then the parameter is set by using an anonymous function, see "help ode45" and Answers: Use anonbymous functions for passing parameters.