Error in using 'frestimate' for MIMO system

2020 6 26
1 个回答
更新时间：2020 6 26
4 次查看（30 天）

0 个投票

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

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

回答（1 个）

Walter Roberson
Walter Roberson 2020-6-26

0 个投票

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
Siddhesh Rane 2020-6-26
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
Walter Roberson 2020-6-26
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
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
Walter Roberson 2020-6-26
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')

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 General Applications 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by