How to set a constraint to limit rate of change of output in s function

I have a simulink model where I have used an s-function block and based on the value of Preq and SOC the s-function provides an output based the algorithm shown below. I want the rate of change of output y(1) to be limited to 50. Please let me know how to do it. I am using fmincon optimization to obtain the output. I think a non linear constraint must be used to do this where the difference between the current output value and the previous output value must be 50. but I do not know how to call the current output value and previous output value while using an s function

 function [sys,x0,str,ts,simStateCompliance] = trail1(t,x,u,flag)

switch flag

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

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

  
  % Update %
  
  case 2
    sys=mdlUpdate(t,x,u);

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

  
  % GetTimeOfNextVarHit %
  
  case 4
    sys=mdlGetTimeOfNextVarHit(t,x,u);

 
  % Terminate %
 
  case 9
    sys=mdlTerminate(t,x,u);

  % Unexpected flags %
  
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));

end

% end sfuntmpl

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


sizes = simsizes;

sizes.NumContStates  = 0;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 3;  %% the three outputs are the power delivered by fuel cell, battery and the equivalent factor  
sizes.NumInputs      = 2;  %% Inputs are Power required and SOC
sizes.DirFeedthrough = 1;  %% because output is a function of SOC 
sizes.NumSampleTimes = 1;  %% at least one sample time is needed

sys = simsizes(sizes);

%
% initialize the initial conditions
%
x0  = [];

%
% str is always an empty matrix
%
str = [];

%
% initialize the array of sample times
%
ts  = [-1 0];

% Specify the block simStateCompliance. The allowed values are:
%    'UnknownSimState', < The default setting; warn and assume DefaultSimState
%    'DefaultSimState', < Same sim state as a built-in block
%    'HasNoSimState',   < No sim state
%    'DisallowSimState' < Error out when saving or restoring the model sim state
simStateCompliance = 'UnknownSimState';

% end mdlInitializeSizes

%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u) %#ok<*INUSD>

sys = [];

% end mdlDerivatives

%
%=============================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=============================================================================
%
function sys=mdlUpdate(t,x,u)

sys = [];

% end mdlUpdate

%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys=mdlOutputs(t,x,u)
%constants initialization
Pbatt_char= 7000; 
SOC_min=0.4;
SOC_max=0.90;
Pfc_min=5000;
Pfc_max=7000;
Pbatt_max=15000;

%define Matrix Aeq

Aeq=[0 1 0;1 0 1];

%define Matrix beq
mu = 0.6;

beq=[(1-2*mu*((u(2)-0.5*(SOC_max-SOC_min))/(SOC_max-SOC_min))); u(1)];

%define boundary conditions
lb=[Pfc_min, 0, -Pbatt_char];
ub=[Pfc_max, 100, Pbatt_max];

    
   
%defining initial conditions

x0 = [7000 0 6000];
options = optimoptions('fmincon','Algorithm','active-set','Display','off','MaxFunctionEvaluations',1000,'MaxIterations',100);
[y,fval] = fmincon('OF_ECMS',x0,[],[],Aeq,beq,lb,ub,[],options); %#ok<*ASGLU>

Pfc=y(1); Pbatt=y(3); alpha=y(2); 
    
    


sys = [Pfc Pbatt alpha];

% end mdlOutputs

%
%=============================================================================
% mdlGetTimeOfNextVarHit
% Return the time of the next hit for this block.  Note that the result is
% absolute time.  Note that this function is only used when you specify a
% variable discrete-time sample time [-2 0] in the sample time array in
% mdlInitializeSizes.
%=============================================================================
%
function sys=mdlGetTimeOfNextVarHit(t,x,u)
sampleTime = 1;    %  Example, set the next hit to be one second later.
sys = t + sampleTime;

% end mdlGetTimeOfNextVarHit

%
%=============================================================================
mdlTerminate
Perform any end of simulation tasks.
% =============================================================================

function sys=mdlTerminate(t,x,u)
    
sys = [];
  

  

% end mdlTerminate

Objective function code is

function f = OF_ECMS(x)

f = (x(1)+x(2)*x(3));

[Mathieu NOE]

hello

why not using the aleady available simulink block ?

Limit rate of change of signal - Simulink - MathWorks France

[aditya deepak]

Thank you for your reply. There is a requirement that the sum of outputs y(1) and y(3) should equal to demand and by using the rate limiter that requirement would not be satisfied. I am using fmincon optimization to obtain the output. I think a non linear constraint must be used to do this where the difference between the current output value and the previous output value must be 50. but I do not know how to call the current output value and previous output value while using an s function.

Answer 1

[Stephen Eltinge]

Hi aditya,

To write an S-Function block that recalls its previous output, please look into adding a discrete state to the S-Function. See the "unit delay" Level-2 MATLAB S-Function example for a simple starting point. You should then be able to use the state's data to provide an additional constraint for your optimization function.

[aditya deepak]

Thank you for your reply

Should I use a Level-2 S-function or is it possible with a Level 1?

[Stephen Eltinge]

I see that the code you posted is for a Level-1 S-Function. Yes, a Level-1 S-Function has the capabilities that you need (see this page for more information). When starting from scratch, though, we generally recommend using Level-2 S-Functions, since they support a fuller range of Simulink features.

[aditya deepak]

Is there a sample code where discrete state in a level1 s function is used?

[Stephen Eltinge]

I don't think we maintain examples for Level-1 S-Functions anymore. You could start with a Level-2 S-Function and use the table here to translate back into Level-1 code. If you need more help, feel free to message MathWorks support!

[aditya deepak]

I converted the level 1 s function to level 2 s function and I am getting the following error

Undefined function 'OutPort' for input arguments of type 'Simulink.MSFcnRunTimeBlock'.

function ECMSsf2(block)

setup(block);

%endfunction

function setup(block)

%% Register number of input and output ports

block.NumInputPorts = 2;

block.NumOutputPorts = 3;

%% Setup functional port properties to dynamically

%% inherited.

block.SetPreCompInpPortInfoToDynamic;

block.SetPreCompOutPortInfoToDynamic;

block.InputPort(1).DirectFeedthrough = false;

%% Set block sample time to inherited

block.SampleTimes = [-1 0];

%% Set the block simStateCompliance to default (i.e., same as a built-in block)

block.SimStateCompliance = 'DefaultSimState';

%% Register methods

block.RegBlockMethod('PostPropagationSetup',@DoPostPropSetup);

block.RegBlockMethod('InitializeConditions',@InitConditions);

block.RegBlockMethod('SetInputPortSamplingMode', @SetInpPortFrameData);

block.RegBlockMethod('Outputs', @Output);

block.RegBlockMethod('Update', @Update);

function DoPostPropSetup(block)

%% Setup Dwork

block.NumDworks = 1;

block.Dwork(1).Name = 'x0';

block.Dwork(1).Dimensions = 1;

block.Dwork(1).DatatypeID = 0;

block.Dwork(1).Complexity = 'Real';

block.Dwork(1).UsedAsDiscState = true;

end

function InitConditions(block)

% Initialize Dwork

block.Dwork(1).Data = 0;

end

function SetInpPortFrameData(block, idx, fd)

block.InputPort(idx).SamplingMode = fd;

block.InputPort(idx).SamplingMode = fd;

block.OutputPort(1).SamplingMode = fd;

block.OutputPort(2).SamplingMode = fd;

block.OutputPort(3).SamplingMode = fd;

end

function Output(block)

%constants initialization

Pbatt_char=5000;

SOC_min=0.2;

SOC_max=0.90;

Pfc_min=1000;

Pfc_max=10000;

Pbatt_max=17000;

%define Matrix Aeq

Aeq=[0 1 0;1 0 1];

%define Matrix beq

mu = 0.6;

beq=[(1-2*mu*((block.InputPort(2).Data-0.5*(SOC_max+SOC_min))/(SOC_max+SOC_min))); block.InputPort(1).Data];

%define boundary conditions

lb=[Pfc_min, 0, -Pbatt_char];

ub=[Pfc_max, 100, Pbatt_max];

%defining initial conditions

x0 = [0 0 0];

options = optimoptions('fmincon','Algorithm','active-set','Display','off','MaxFunctionEvaluations',1000,'MaxIterations',100);

nonlcon=@(y) Nonlinearequations(y,block);

[y,fval] = fmincon(@OF_ECMS,x0,[],[],Aeq,beq,lb,ub,nonlcon,options);

function f = OF_ECMS(y) %<----fix

f = (y(1)+y(2)*y(3));

end

function [c, ceq] = Nonlinearequations(y,block)%<----fix

c = [];

ceq = [];

PFCold = block.Dwork(1).Data;

c(1) = PFCold - y(1) - 20; %<----change to -20

c(2) = y(1) - PFCold - 20;

end

Pfc=y(1); Pbatt=y(3); alpha=y(2);

block.OutPort(1).Data = Pfc;

block.OutPort(2).Data = Pbatt;

block.OutPort(3).Data = alpha;

end

function Update(block)

block.Dwork(1).Data = block.Outport(1).Data;

end

end

end

[Stephen Eltinge]

Hi aditya,

Please try replacing "OutPort" with "OutputPort" in its four occurrences near the end of your code.