Error in using 'frestimate' for MIMO system

Question

Siddhesh Rane 2020-6-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/555142-error-in-using-frestimate-for-mimo-system

评论： Walter Roberson 2020-6-26

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2020-6-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/555142-error-in-using-frestimate-for-mimo-system#answer_457423

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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Siddhesh Rane 2020-6-26

在 MATLAB Online 中打开

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  

Siddhesh Rane 2020-6-26

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 2020-6-26

在 MATLAB Online 中打开

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')

请先登录，再进行评论。