level-2 matlab s-function error
Show older comments
Hi,
please, I'd like to know what's wrong with the following code. When I run it, I get the following error message: "Derivative of state '1' in block 'model_4WID/Level-2 MATLAB S-Function3' at time 0.0 is not finite. The simulation will be stopped. There may be a singularity in the solution. If not, try reducing the step size (either by reducing the fixed step size or by tightening the error tolerances)"
I really need your help.
Thank you.
function dynamic(block)
setup(block);
function setup(block)
% number of ports.
block.NumInputPorts = 1;
block.NumOutputPorts = 1;
% port properties.
block.SetPreCompInpPortInfoToDynamic;
block.SetPreCompOutPortInfoToDynamic;
%input port properties.
block.InputPort(1).DatatypeID = 0; % double
block.InputPort(1).Complexity = 'Real';
block.InputPort(1).Dimensions = 20;
block.InputPort(1).DirectFeedthrough = false;
%output port properties.
block.OutputPort(1).DatatypeID = 0; % double
block.OutputPort(1).Complexity = 'Real';
block.OutputPort(1).Dimensions = 18;
% continuous states.
block.NumContStates = 5;
block.SampleTimes = [0 0];
% Specify the block simStateCompliance. The allowed values are:
% 'DefaultSimState', < Same SimState as a built-in block
block.RegBlockMethod('InitializeConditions', @InitializeConditions);
% block.RegBlockMethod('Start', @Start);
block.RegBlockMethod('Outputs', @Outputs);
block.RegBlockMethod('Derivatives', @Derivatives);
%endfunction
function InitializeConditions(block)
block.ContStates.Data(1) = 20;
block.ContStates.Data(2)= 0;
block.ContStates.Data(3)= 0;
block.ContStates.Data(4)= 0;
block.ContStates.Data(5)= 0;
%endfunction
function Outputs(block)
de=zeros(1,5); u=zeros(1,20);
for i=1:5
de(i)=block.ContStates.Data(i);
end
for i=1:22
u(i)=block.InputPort(1).data(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% non-linear dugoff tyre model %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cz=0.001; %vertical deflection rate of the tyre
Cs=50000;
epsilon=0.015;
I=2.1;
Iz=1627;
R=0.35;
Ca=30000;
m=1298.9;
lf=1;
lr=1.454;
br=1.436;
bf=br;
mu=0.9;
g=9.81;
cons=m/(lr+lf);
h=0.5;
r=de(3); phi=atan2(de(2),de(1));
C=[0 m*r 0 0 0;-m*r 0 0 0 0;0 0 0 0 0;0 0 0 cos(phi) -sin(phi);0 0 0 sin(phi) cos(phi)];
M=[m 0 0 0 0;0 m 0 0 0;0 0 Iz 0 0;0 0 0 1 0;0 0 0 0 1];
u1=(de(2)+lf*de(3))/(de(1)+0.5*bf*de(3));
u2=(de(2)+lf*de(3))/(de(1)-0.5*bf*de(3));
u3=(de(2)-lr*de(3))/(de(1)+0.5*br*de(3)); %% les ui représentent les rapports permettant de calculer les angles de dérives de chaque roue%%
u4=(de(2)-lr*de(3))/(de(1)-0.5*br*de(3));
U=[u1 u2 u3 u4];
%------vertical load/charge verticale:----------
dde=(-de'*C)/M; %dérivée de de!
Fz1= cons*(0.5*g*lr-0.5*dde(1)*h-lr*h*dde(2)/bf);
Fz2= cons*(0.5*g*lr-0.5*dde(1)*h+lr*h*dde(2)/bf);
Fz3= cons*(0.5*g*lr+0.5*dde(1)*h-lr*h*dde(2)/bf);
Fz4= cons*(0.5*g*lr+0.5*dde(1)*h+lr*h*dde(2)/bf);
Fz=[Fz1 Fz2 Fz3 Fz4]';
%%% velocity component in the wheel plane : is the longitunal velocity
v1=(de(1)+0.5*bf*de(3))*cos(u(1))+(de(2)+lf*de(3))*sin(u(1));
v2=(de(1)-0.5*bf*de(3))*cos(u(2))+(de(2)+lf*de(3))*sin(u(2));
v3=(de(1)+0.5*br*de(3))*cos(u(3))+(de(2)-lr*de(3))*sin(u(3));
v4=(de(1)-0.5*br*de(3))*cos(u(4))+(de(2)-lr*de(3))*sin(u(4));
V=[v1 v2 v3 v4]';
lamda=zeros(4,1);f=length(lamda);
omega=zeros(4,1);
alph=zeros(4,1);
Re=zeros(4,1);
S=zeros(4,1);
dz=zeros(4,1);
Fs=zeros(4,1);
Ft=zeros(4,1);
for i=1:4
dz(i)=-Cz*Fz(i)+ 0.33*R;
Re(i)=R-dz(i)./3;
omega(i)=(-R*u(i+16)+u(i+4))/I;
S(i)=1+omega(i)*Re(i)./V(i);
alph(i)=atan(U(i))-u(i);
lamda(i)= mu*Fz(i)*(1-S(i))*(1-epsilon*V(i)*sqrt((S(i))^2 + (tan(alph(i)))^2))/(2*sqrt((Cs^2)*(S(i))^2 + (Ca^2)*(tan(alph(i))^2)));
if lamda(i)<1
f(i)=lamda(i)*(2-lamda(i));
Fs(i)=Ca*f(i)*tan(alph(i))/(1-S(i));
Ft(i)=Cs*f(i)*S(i)/(1-S(i));
elseif lamda(i)>1
f(i)=1;
Fs(i)=Ca*f(i)*tan(alph(i))/(1-S(i));
Ft(i)=Cs*f(i)*S(i)/(1-S(i));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Sorties %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
block.OutputPort(1).Data(1)=de(1);
block.OutputPort(1).Data(2)=de(2);
block.OutputPort(1).Data(3)= phi;
block.OutputPort(1).Data(4)=de(4);
block.OutputPort(1).Data(5)=de(5);
block.OutputPort(1).Data(6)=Ft(1);
block.OutputPort(1).Data(7)=Ft(2);
block.OutputPort(1).Data(8)=Ft(3);
block.OutputPort(1).Data(9)=Ft(4);
block.OutputPort(1).Data(10)=Fs(1);
block.OutputPort(1).Data(11)=Fs(2);
block.OutputPort(1).Data(12)=Fs(3);
block.OutputPort(1).Data(13)=Fs(4);
block.OutputPort(1).Data(14)= r;
block.OutputPort(1).Data(15)=Fz(1);
block.OutputPort(1).Data(16)=Fz(2);
block.OutputPort(1).Data(17)=Fz(3);
block.OutputPort(1).Data(18)=Fz(4);
%endfunction
function Derivatives(block)
m=1298.9; Iz=1627;lf=1; lr=1.454; br=1.436; bf=br;
de=zeros(1,5); u=zeros(1,22);
for i=1:5
de(i)=block.ContStates.Data(i);
end
for i=1:20
u(i)=block.InputPort(1).Data(i);
end
r=de(3); phi=atan2(de(2),de(1));
M=[m 0 0 0 0;0 m 0 0 0;0 0 Iz 0 0;0 0 0 1 0;0 0 0 0 1]; %matrice d'inertie
C=[0 m*r 0 0 0;-m*r 0 0 0 0;0 0 0 0 0;0 0 0 cos(phi) -sin(phi);0 0 0 sin(phi) cos(phi)];
F=[u(9)*cos(u(1))-u(13)*sin(u(1))+u(10)*cos(u(2))-u(14)*sin(u(2))+u(11)*cos(u(3))-u(15)*sin(u(3))+u(12)*cos(u(4))-u(16)*sin(u(4));u(13)*cos(u(1))+u(9)*sin(u(1))+u(14)*cos(u(2))+u(10)*sin(u(2))+u(15)*cos(u(3))+u(11)*sin(u(3))+u(16)*cos(u(4))+u(12)*sin(u(4));u(9)*(lf*sin(u(1))+0.5*bf*cos(u(1)))+u(10)*(lr*sin(u(2))-0.5*bf*cos(u(2)))+u(11)*(-lr*sin(u(3))+0.5*br*cos(u(3)))+u(12)*(-lr*sin(u(4))-0.5*br*cos(u(4)))+u(13)*(lf*sin(u(1))-0.5*bf*cos(u(1)))+u(14)*(lf*sin(u(2))+0.5*bf*cos(u(2)))+u(15)*(-lr*sin(u(3))-0.5*br*cos(u(3)))+u(16)*(-lr*sin(u(4))+0.5*br*cos(u(4)));0;0];
X=(F'-C*de)/M;
block.Derivatives.Data(1) = X(1);
block.Derivatives.Data(2) = X(2);
block.Derivatives.Data(3) = X(3);
block.Derivatives.Data(4) = X(4);
block.Derivatives.Data(5) = X(5);
%endfunction
Answers (1)
Fangjun Jiang
on 24 Feb 2020
This line is invalid. Dimension mismatch
X=(F'-C*de)/M
1 Comment
Valéry Ebogo
on 6 Mar 2020
Categories
Find more on Startup and Shutdown in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!