level-2 matlab s-function error

5 views (last 30 days)
Valéry Ebogo
Valéry Ebogo on 15 Feb 2020
Commented: Valéry Ebogo on 6 Mar 2020
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
Fangjun Jiang on 24 Feb 2020
This line is invalid. Dimension mismatch
X=(F'-C*de)/M

Products


Release

R2015a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!