Error in using 'frestimate' for MIMO system

I am using frestimate for finding frequency response function of a non-linear system which has multiple inputs and multiple outputs,

Picture above shows the simulink model, I have selected 1 input (out of 2) and 1 output (output), I get following error after using

 [sysest,simout] = frestimate(mdl,getlinio(mdl),in);

Error using frestimate>LocalRunSimulation (line 1165)
State derivatives returned by S-function 'System_Model' in 'Non_Linear_Model_FRF/Validation
Setup/S-Function' during flag=1 call must be a real vector of length 13

Error in frestimate (line 255)
            [simout{ctexp},err] = LocalRunSimulation(parammgr,SimulationPackage,ctexp); -
            Show complete stack trace

Answer 1

[Walter Roberson]

Your S function Validation Setup has 13 outputs. The model derivatives for it must return 13 values, one for each of those 13 outputs.

It looks to me as if your S function is probably a Level 1 S function that is being passed u and flag and is expected to detect the various flag states to decide what it is doing.

[Siddhesh Rane]

Hi Walter, Thanks for the response. I didn't quite understand it. Do I have to mark all 13 outputs in the simulink model as "Output Measurement" ?

I want to estimate frequency response function for 1 output (out of 13 )and 1 input (out of 2). Below is the code for my s-function, as you pointed out it's level 1.

function [sys,x0,str,ts] = System_Model(t,x,u,flag)
%% S- Function for model of validation setup
%% Last Updated : 24 June 2020
%% last Updated : Modified for variable radii, 2 additional states added for the radii of 

%% Inputs u; 
%u(1)= tau1
%u(2)= tau4

%% Outputs y
%y(1)=x(1)
%y(2)=x(2)


switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,
    [sys,x0,str,ts]=mdlInitializeSizes;

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,
    sys=mdlDerivatives(t,x,u);

  %%%%%%%%%%%
  % Outputs %
  %%%%%%%%%%%
  case 3,
    sys=mdlOutputs(t,x,u);

  %%%%%%%%%%%%%%%%%%%
  % Unhandled flags %
  %%%%%%%%%%%%%%%%%%%
  case { 2, 4, 9 },
    sys = [];

  %%%%%%%%%%%%%%%%%%%%
  % Unexpected flags %
  %%%%%%%%%%%%%%%%%%%%
  otherwise
    error(['Unhandled flag = ',num2str(flag)]);

end
% end wpfun1

%
%=============================================================================
%% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts]=mdlInitializeSizes

sizes = simsizes;
sizes.NumContStates  = 13;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 13;
sizes.NumInputs      = 2;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;

sys = simsizes(sizes);
x0  = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0.250; 0.0375];
str = [];
ts  = [0 0];

% end mdlInitializeSizes
%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u)
        
%% Model Parameters

% Initial roll radius of unwinder in m
% R0 = 0.3;

% Initial roll radius of Rewinder in m
% R04 = 0;

% Coeefficient of viscous friction for bearings (N-m-s/rad)
b = 0;

% Moment of inertia of the spool (kg-m2)
J_spool = 0.05; 

% Density of the tape (kg/m3)
tape_density = 1.46e3;

% Width of the tape (m)
tape_width = 18e-3;

% Thickmess of the tape (m)
tape_thickness = 0.19e-3;

% Elastic modulus of the paper 
E = 3e9;

% Outer radius of the spool/ Inner radius of the tape (m)
Rc = 37.5e-3;

% Radius of idle roller 1 (m)
R2 = 25e-3;

% Radius of idle roller 2 (m)
R3 = 25e-3;

% Length of the span 1 (unwinding roller - idle roller 1)
% L1 = sqrt(d1^2 - (R1-R2)^2)
% d1 = distance between centers of un-winder spool and idler roller 1 = 0. 59471 m
d1 = 0.59471;


% Length of the span 2 (idle roller 1 - idle roller 2)
L2 = 1.5;

% Length of the span 3 (idle roller 2 - rewinding roller)
% L3 = sqrt(d3^2 - (R4-R3)^2)
% d3 = distance between centers of re-winder spool and idler roller 2 = 0. 59471 m
d3 = 0.59471;

%Moment of inertia of idle roller 1 (kg-m2)
J2 = 6e-4; 

%Moment of inertia of idle roller 2 (kg-m2)
J3 = 6e-4; 

%% Unwinder radius estimation
% R1 = R0 -((y(x1)/2*pi)*tape_thickness);

%% Unwinder inertia calculations
% J1 = Jspool + ((pi/2)*(tape_width*tape_density)*((R1)^4-(Rc)^4));

%% Rewinder radius estimation
% R4 = R04 + ((y(x7)/2*pi)*tape_thickness);

%% Unwinder inertia calculations
% J1 = Jspool + ((pi/2)*(tape_width*tape_density)*((R1)^4-(Rc)^4));

%% Inputs
tau1 =u(1);			%oC
tau4= u(2);			%oC

%% System Model

% Unwinder force balance
xdot(1) = x(2);
xdot(2) = (x(12)*x(9) - tau1 + ...
    (tape_density*tape_width*tape_thickness*(x(12)^3)*(x(2)^2)))...
    /(J_spool + ((pi/2)*(tape_width*tape_density)*((x(12))^4-(Rc)^4)));

% Rewinder force balance
xdot(7) = x(8);
xdot(8) = (-x(13)*x(11) + tau4 - ...
    (tape_density*tape_width*tape_thickness*(x(13)^3)*(x(8)^2)))...
    /(J_spool + ((pi/2)*(tape_width*tape_density)*((x(13))^4-(Rc)^4)));


% Idle roller 1 force balance
xdot(3) = x(4);
xdot(4) = (R2*(x(10)-x(9)))/J2;


% Idle roller 2 force balance
xdot(5) = x(6);
xdot(6) = (R3*(x(11)-x(10)))/J3;

% Mass balance for span 1 (unwinding roller to idle roller 1)
xdot(9) = (((E*tape_width*tape_thickness)*((R2*x(4))-(x(12)*x(2)))) + (-x(9)*R2*x(4)))/(sqrt(d1^2 - (x(12)-R2)^2));

% Mass balance for span 2 (idle roller 1 to idle roller 2)
xdot(10) = ((E*tape_width*tape_thickness)*(R3*x(6)-R2*x(4)) + (R2*x(4)*x(9)-R3*x(6)*x(10)))/L2;  

% Mass balance for span 3 (idle roller 3 to rewinding roller)
xdot(11) = ((E*tape_width*tape_thickness)*(x(13)*x(8)-R3*x(6)) + (R3*x(6)*x(10)-x(13)*x(11)*x(8)))/sqrt(d3^2 - (x(13)-R3)^2); 

% Un-winder roll radius change
xdot(12) = -(tape_thickness* x(2))/(2*pi);

% Un-winder roll radius change
xdot(13) = (tape_thickness* x(8))/(2*pi);


sys = [xdot(1); xdot(2); xdot(3); xdot(4); xdot(5); xdot(6); xdot(7); xdot(8); xdot(9); xdot(10); xdot(11); xdot(12); xdot(13)];

% end mdlDerivatives
%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys=mdlOutputs(t,x,u)

sys = x;


% end mdlOutputs  

[Walter Roberson]

In calculationg xdot(9) and xdot(11) you have some sqrt() of subtractions. Those subtractions can come out negative, leading to sqrt() of negative numbers, leading to non-real results.

[Siddhesh Rane]

I changed it to real-positive values (0.6), which did not solve the problem. I am not able to understand from error message if the problem is with non-real values or length of the output vector. Because I am selecting only 1 input and 1 output.

[Walter Roberson]

After the assignment to sys, add in lines:

assert(length(sys) == 13, 'derivatives are wrong length')
assert(all(isfinite(sys)), 'derivatives are not all finite')
assert(all(imag(sys) == 0), 'derivatives are not all real')