During flag=1 call must be a real vector of length 10 error

Hi all

I have a nonlinear matlab code which is

function[sys,x0,str,ts]=nonlinearxrae1_height(t,x,ui,flag,x0,s,ar,zf1,zf2,b,et,lf1,lf2,c,ha,mass,ix,iy,iz,ixz,clow,clo,clq,cdow,cmo,cmq,claw,cla,cln,cma,cmn,clad,k,cmad,yvb,lx,yzb,lzb,nzb,ro,g,pi)
%An M-file S-function for defining a system of continuous state equations:
%XRAE1 6DOF nonlinear simulation
%Dispatch the flag.

switch flag,
    case 0
        [sys,x0,str,ts]=mdlInitializeSizes(x0)
    case 1
        sys=mdlDerivatives(t,x,ui,s,ar,zf1,zf2,b,et,lf1,lf2,c,ha,mass,ix,iy,iz,ixz,clow,clo,clq,cdow,cmo,cmq,claw,cla,cln,cma,cmn,clad,k,cmad,yvb,lx,yzb,lzb,nzb,ro,g,pi)
    case 3
        sys=mdlOutputs(t,x); %Calculate outputs
    case{2,4,9} %Unused flags
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]); %Error handling
end

%End of XRAE1 func.
%=========================================================================
%mdlInitializesSizes
%Return the sizes, inital conditions and sample times for the S-function.
function [sys,x0,str,ts]=mdlInitializeSizes(x0)

%call simsizes for a sizes structure, fill it in and convert it
%to a sizes array.
sizes=simsizes;

sizes.NumContStates=10;
sizes.NumDiscStates=0;
sizes.NumOutputs=10;
sizes.NumInputs=4;
sizes.DirFeedthrough=0; %matrix Dis nonempty
sizes.NumSampleTimes=1;

sys=simsizes(sizes);

%initialize the inital conditions
x0=x0;
%str is an empty matrix.
str=[];
%Initialize the array of sample times; in this example the sample
%time is continuous, so set ts to 0 and its offset to 0.
ts=[0 0];
%End of mdlInitializeSizes.

%========================================================================
%mdlDerivatives
%Return the derivatives for the continuous states.
%========================================================================

function sys=mdlDerivatives(t,x,ui,s,ar,zf1,zf2,b,et,lf1,lf2,c,ha,mass,ix,iy,iz,ixz,clow,clo,clq,cdow,cmo,cmq,claw,cla,cln,cma,cmn,clad,k,cmad,yvb,lx,yzb,lzb,nzb,ro,g,pi)
%angle of attack
A=atan(x(3)/x(1));
%total velocity
vt=sqrt(x(1)^2+x(2)^2+x(3)^2);
%dynamic pressure
qp=0.5*ro*(vt^2);
%part of total lift coeficient
cl1=clo+cla*A+clq*((x(5)*c)/(2*vt))+cln*ui(1);
%wing body lift coefficent
clw=clow+claw*A;
%total drag coefficient
cd=cdow+k*(clw^2);
%thrust
T=26.7154*ui(4)-0.0026*(vt^2);

MAS=[1+qp*s*x(3)*clad*c*sin(A)/(2*mass*vt*(x(1)^2+x(3)^2)),-qp*s*x(1)*clad*c*sin(A)/(2*mass*vt*(x(1)^2+x(3)^2));-qp*s*x(3)*clad*c*cos(A)/(2*mass*vt*(x(1)^2+x(3)^2)),1+qp*s*x(1)*clad*c*cos(A)/(2*mass*vt*(x(1)^2+x(3)^2))];
Forces=[x(6)*x(2)-x(5)*x(3)-g*sin(x(8))+qp*s*(cl1*sin(A)-cd*cos(A))/mass+T/mass;x(5)*x(1)-x(4)*x(2)+g*cos(x(8))*cos(x(7))-qp*s*(cl1*cos(A)+cd*sin(A))/mass];
%[sys(1),sys(3)]=inv(MAS)*Forces;

x13d=MAS\Forces;
sys(1)=x13d(1);
sys(3)=x13d(2);
% coeficient needed for longitudinal equation
%b1=x(6)*x(2)-x(5)*x(3)-g*sin(x(8))+qp*s*(cl1*sin(A)-cd*cos(A))/mass+T/mass;
%b2=x(5)*x(1)-x(4)*x(2)+g*cos(x(8))*cos(x(7))-qp*s*(cl1*cos(A)+cd*sin(A))/mass;
%d=qp*s*clad*c/(2*mass*vt*(x(1)^2+x(3)^2));
%d1=1+d*(x(1)*cos(A)+x(3)*sin(A));
%d2=d*(b1*cos(A)+b2*sin(A));
%velocity rate along xb axis
%sys(1)=(b1+x(1)*d2)/d1;
%veloctiy rate along zb axis
%sys(3)=(b2+x(3)*d2)/d1;
%angle of attack rate
ad=(sys(3)*x(1)-sys(1)*x(3))/(x(1)^2+x(3)^2);
%pitching moment coefficient
cm=cmo+cma*A+((cmad*ad*c)/(2*vt))+((cmq*x(5)*c)/(2*vt))+cmn*ui(1);
%lift coefficient
cl=cl1+((clad*ad*c)/(2*vt));
%pitching angular acceleration about yb axis
sys(5)=((iz-ix)*x(4)*x(6)-(x(4)-x(6))*ixz+qp*s*c*cm-T*et-ha*qp*s*(cl*sin(A)-cd*cos(A)))/iy;
%lateral stability dervatives
%stability axes
nv=-0.0363+0.1969*(lf1*cos(A)+zf1*sin(A))/b;
lv=-0.0119-0.0074*clw-0.1969*(zf1*cos(A)-lf1*sin(A))/b;
x1=(zf2-(zf2*cos(A)-lf2*sin(A)))/b;
%sidewash component of ypf
X=[0.0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25];
Y=[0.0 0.036 0.075 0.114 0.156 0.2 0.25 0.295 0.345 0.4 0.45];

y1=abs(x1);
ysd=interp1(X,Y,y1);
yps=x1*ysd/y1;
ypf=-0.3133*((zf2*cos(A)-lf2*sin(A))/b-0.18-yps);
yp=0.078*cl+ypf;
cdvd=2*clw*(k-1/(pi*ar))*claw*pi/180.0;
np=-0.034*cl+1.23*cdvd-ypf*(lf2*cos(A)+zf2*sin(A))/b;
lp=-0.2457+ypf*(zf2*cos(A)-lf2*sin(A))/b;
yrf=0.2164*(lf1*cos(A)+zf1*sin(A))/b;
yr=-0.0109+yrf;
nr=-0.0022-0.1621*cdow-0.009*clw^2-yrf*(lf1*cos(A)+zf1*sin(A))/b;
lr=-0.00189+0.1243*clw+yrf*(zf1*cos(A)-lf1*sin(A))/b;
nx=0.0195*clw;
%transformation to body reference axes
ypb=yp*cos(A)-yr*sin(A);
yrb=yr*cos(A)+yp*sin(A);
lvb=lv*cos(A)-nv*sin(A);
lpb=lp*(cos(A))^2-(lr+np)*sin(A)*cos(A)+nr*(sin(A))^2;
lrb=lr*(cos(A))^2-(nr-lp)*sin(A)*cos(A)-np*(sin(A))^2;
nvb=nv*cos(A)+lv*sin(A);
npb=np*(cos(A))^2-(nr-lp)*sin(A)*cos(A)-lr*(sin(A))^2;
nrb=nr*(cos(A))^2+(lr+np)*sin(A)*cos(A)+lp*(sin(A))^2;
lxb=lx*cos(A)-nx*sin(A);
nxb=nx*cos(A)+lx*sin(A);
%aerodynamic force along yb axis
y=qp*s*(yvb*x(2)+b*ypb*x(4)+b*yrb*x(6))/vt+qp*s*yzb*ui(3);
%aerodynamic moment about xb axes
l=qp*s*b*(lvb*x(2)+b*lpb*x(4)+b*lrb*x(6))/vt+qp*s*b*(lzb*ui(3)+lxb*ui(2));
%aerodynamic moment about zb axis
n=qp*s*b*(nvb*x(2)+b*npb*x(4)+b*nrb*x(6))/vt+qp*s*b*(nzb*ui(3)+nxb*ui(2));
%velocity rate along yb axes
sys(2)=x(4)*x(3)-x(6)*x(1)+g*cos(x(8))*sin(x(7))+y/mass;
Inertia=[ix -ixz; -ixz iz];
Ml=l-(x(5)*x(6)*(iz-iy))+(x(4)*x(5)*ixz);
Mn=n-(x(4)*x(6)*(ix-iz))-((x(4))^2-((x(6))^2)*ixz);
%[sys(4),sys(6)]=inv(Inertia)*[Ml;Mn];
x46d=Inertia\[Ml;Mn];
sys(4)=x46d(1);
sys(6)=x46d(2);

%rolling angular acceleration about zb axes
%sys(4)=((iy-iz)*x(5)*x(6)+l)/ix;
%yawing angular acceleration about zb axes
%sys(6)=((ix-iy)*x(4)*x(5)+n)/iz;
%euler angle rates
sys(7)=x(4)+x(5)*tan(x(8))*sin(x(7))+x(6)*tan(x(8))*cos(x(7));
sys(8)=x(5)*cos(x(7))-x(6)*sin(x(7));
sys(9)=(x(6)*cos(x(7))+x(5)*sin(x(7)))/cos(x(8));
sys(10)=x(1)*sin(x(8))-x(2)*sin(x(7))*cos(x(8))-x(3)*cos(x(7))*cos(x(8));
%end of "dynamic derivative"
%end of program
%End of mdlDerivatives.

%=========================================================================
%mdlOutputs
%Return the block outputs.
%=========================================================================


function sys=mdlOutputs(t,x)
sys=x;

    

This is the code I'm trying to implement in Simulink s file with the s file being the name of the function above.

When I run this I get the error and I don't know why. Any suggestions would be very helpful.

[Maaz Madha]

This is the whole error

[Walter Roberson]

Notice the NaN in the outputs.

As outside observers, we have no way of knowing that your mass is non-zero or that the results of the \ will definitely be finite.

[Maaz Madha]

The parameters for this code are in another code which is run first

%RPV GEOMETRY
s=.9307;
ar=7.48;
zf1=0.0887;
zf2=.1132;
b=2.638;
et=.16;
lf1=1.0951;
lf2=1.1019;
c=.353;
ha=.045;

%RPV MASS & INERTIA CHARACTERISTICS
mass=18.5;
ix=5.0;
iy=2.10;
iz=5.80;
ixz=0.17;


%RPV STABILITY & CONTROL DERIVATIVES
clow=0.79;
clo=0.79;
claw=4.43;
cla=4.97;
clad=2.86;
clq=5.25;
cln=0.48;
cdow=0.02245;
k=0.0520;
cmo=-0.11;
cma=-1.24;
cmad=-9.56;
cmq=-21.27;
cmn=-1.62;
yvb=-0.3054;
lx=-0.2291;
yzb=0.1393;
lzb=0.00468;
nzb=-0.0579;
ro=1.225;
g=9.80665;

% Additional parameters for Model
%alpha = -0.0867;
%vt0=30;
%U = vt0*cos(alpha);
%W = vt0*sin(alpha);

so the mass wouldn't be 0 but is 18.5 and thus it shouldn't be a non zero quantity

[Maaz Madha]

It also worked previously but only now its not working

[Walter Roberson]

If you can attach your .slx then someone might be willing to test.

[Maaz Madha]

I have attached the slx cw

The K_xrae is

or

K_xrae=zeros(4,10);
K_xrae(1,1)=-0.2488;
K_xrae(1,3)=0.3878;
K_xrae(1,5)=-0.4488;
K_xrae(1,8)=-19.9541;
K_xrae(1,10)=-2.3526;
K_xrae(4,1)=15.2227;
K_xrae(4,3)=-9.1902;
K_xrae(4,5)=5.2662;
K_xrae(4,8)=375.5046;
K_xrae(4,10)=29.4992;

Answer 1

[Chetan]

Hello @Maaz Madha,

I Understand that you're working on a MATLAB S-function for a 6DOF nonlinear system. And facing issues while simulating the model.

It seems like you are using level-1 S-functions. You can follow these steps to solve the issues:

Update Fortran S-Function: Level-1 Fortran S-functions are deprecated.

Transition to a Level-2 S-Function using a C gateway as outlined here:

https://www.mathworks.com/help/simulink/sfg/creating-level-2-fortran-s-functions.html

MATLAB S-function Troubleshooting:

• Script Name and Path: Ensure `nonlinearxrae1_height.m` is named correctly and in MATLAB's path.
• Input Connections: Check that all Simulink inputs are connected and dimensioned correctly.
• Initial Conditions: The `x0` vector should correspond to your model's state count (10).
• Variable Initialization: Confirm initialization and passing of all variables and parameters.
• Error Diagnosis: Look at MATLAB command window errors for guidance.
• Debugging: Use breakpoints to debug your S-function script.
• S-Function Registration: Use the correct S-Function block name in Simulink (`nonlinearxrae1_height`)

Refer to the following MathWorks documentation regarding the level 2 S-Functions:

https://www.mathworks.com/help/simulink/sfg/writing-level-2-matlab-s-functions.html

Hope It helps