Problem in a disturbance rejection tuning controller

Hi guys. I want to tune a controller so it can reject the disturbance of an Input. This is my model's code that I'm using in S-function:

function dxdsdp = fermentor_2(t,x,u)

%Output variables
X = x(1);   %Biomass[g/L]
S = x(2);   %Substrate[g/L]
P = x(3);   %Product concentration[g/L]

%Input variables
D = u(1);   %Dilution rate(D=F/V)[1/h]    %F and V is feed rate and volume of the reactor [F] = L/h , [V] = L
S_f = u(2); %Feed concentration of substrate[g/L]


%Constant Parameters
Y_XS = 0.4;  %Mass of new cells/mass of substrate consumed [-]
betta = 0.2; %Constant[1/h]
P_m = 50;    %Maximum product concentration base?[g/L]
K_i = 22;    %Inhibition constant[g/L]
alpha = 2.2; %Constant[g/g]
miu_m = 0.48;%Maximum specific growth rate[1/h]
K_m = 1.2;   %Substrate saturation constant[g/L]


%Equations
miu = ((miu_m*(1-P/P_m))*S)/(K_m+S+S^2/K_i);    %Growth term[1/h]
dX = -D*X+miu*X;                                %Biomass mass balance
dS = D*(S_f-S)-miu*X/Y_XS;                      %Substrate mass balance
dP = -D*P+(alpha*miu+betta)*X;                  %Product mass balance

%Output
dxdsdp = [dX;dS;dP];
end

S_f is my input that is going to change 10% and I want my controller to reject this change. this is my s-function code:

function [sys,x0,str,ts,simStateCompliance] = fermentor_sfcn_D_2(t,x,u,flag,init)
%SFUNTMPL General MATLAB S-Function Template
%   With MATLAB S-functions, you can define you own ordinary differential
%   equations (ODEs), discrete system equations, and/or just about
%   any type of algorithm to be used within a Simulink block diagram.
%
%   The general form of an MATLAB S-function syntax is:
%       [SYS,X0,STR,TS,SIMSTATECOMPLIANCE] = SFUNC(T,X,U,FLAG,P1,...,Pn)
%
%   What is returned by SFUNC at a given point in time, T, depends on the
%   value of the FLAG, the current state vector, X, and the current
%   input vector, U.
%
%   FLAG   RESULT             DESCRIPTION
%   -----  ------             --------------------------------------------
%   0      [SIZES,X0,STR,TS]  Initialization, return system sizes in SYS,
%                             initial state in X0, state ordering strings
%                             in STR, and sample times in TS.
%   1      DX                 Return continuous state derivatives in SYS.
%   2      DS                 Update discrete states SYS = X(n+1)
%   3      Y                  Return outputs in SYS.
%   4      TNEXT              Return next time hit for variable step sample
%                             time in SYS.
%   5                         Reserved for future (root finding).
%   9      []                 Termination, perform any cleanup SYS=[].
%
%
%   The state vectors, X and X0 consists of continuous states followed
%   by discrete states.
%
%   Optional parameters, P1,...,Pn can be provided to the S-function and
%   used during any FLAG operation.
%
%   When SFUNC is called with FLAG = 0, the following information
%   should be returned:
%
%      SYS(1) = Number of continuous states.
%      SYS(2) = Number of discrete states.
%      SYS(3) = Number of outputs.
%      SYS(4) = Number of inputs.
%               Any of the first four elements in SYS can be specified
%               as -1 indicating that they are dynamically sized. The
%               actual length for all other flags will be equal to the
%               length of the input, U.
%      SYS(5) = Reserved for root finding. Must be zero.
%      SYS(6) = Direct feedthrough flag (1=yes, 0=no). The s-function
%               has direct feedthrough if U is used during the FLAG=3
%               call. Setting this to 0 is akin to making a promise that
%               U will not be used during FLAG=3. If you break the promise
%               then unpredictable results will occur.
%      SYS(7) = Number of sample times. This is the number of rows in TS.
%
%
%      X0     = Initial state conditions or [] if no states.
%
%      STR    = State ordering strings which is generally specified as [].
%
%      TS     = An m-by-2 matrix containing the sample time
%               (period, offset) information. Where m = number of sample
%               times. The ordering of the sample times must be:
%
%               TS = [0      0,      : Continuous sample time.
%                     0      1,      : Continuous, but fixed in minor step
%                                      sample time.
%                     PERIOD OFFSET, : Discrete sample time where
%                                      PERIOD > 0 & OFFSET < PERIOD.
%                     -2     0];     : Variable step discrete sample time
%                                      where FLAG=4 is used to get time of
%                                      next hit.
%
%               There can be more than one sample time providing
%               they are ordered such that they are monotonically
%               increasing. Only the needed sample times should be
%               specified in TS. When specifying more than one
%               sample time, you must check for sample hits explicitly by
%               seeing if
%                  abs(round((T-OFFSET)/PERIOD) - (T-OFFSET)/PERIOD)
%               is within a specified tolerance, generally 1e-8. This
%               tolerance is dependent upon your model's sampling times
%               and simulation time.
%
%               You can also specify that the sample time of the S-function
%               is inherited from the driving block. For functions which
%               change during minor steps, this is done by
%               specifying SYS(7) = 1 and TS = [-1 0]. For functions which
%               are held during minor steps, this is done by specifying
%               SYS(7) = 1 and TS = [-1 1].
%
%      SIMSTATECOMPLIANCE = Specifices how to handle this block when saving and
%                           restoring the complete simulation state of the
%                           model. The allowed values are: 'DefaultSimState',
%                           'HasNoSimState' or 'DisallowSimState'. If this value
%                           is not speficified, then the block's compliance with
%                           simState feature is set to 'UknownSimState'.


%   Copyright 1990-2010 The MathWorks, Inc.

%
% The following outlines the general structure of an S-function.
%
switch flag,

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

  %%%%%%%%%%%%%%%
  % 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(init)

%
% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
%
% Note that in this example, the values are hard coded.  This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
%
sizes = simsizes;

sizes.NumContStates  = 3;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 3;
sizes.NumInputs      = 2;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;   % at least one sample time is needed

sys = simsizes(sizes);

%
% initialize the initial conditions
%
x0  = init;

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

%
% initialize the array of sample times
%
ts  = [0 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)

sys = fermentor_2(t,x,u);

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

sys = x;

% 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

and this is the initial conditions for outputs:

[6.00083606242659 4.99790984393358 44.0544669076550]

Also I attached my simulink file too.(I'm using 2021b and I have exported my file to 2018b for people who have older versions of MATLAB, If you guys run into a problem with opening the simulink file please let me know).

My first input(D) is at 0.0389 and my second input(S_f) is going to change at 0 from 20 to 22. Also this the picture of my simulation.

First my model is not converging and simulink can't tune controller properly. What's wrong with my model and my simulation? I haven't work with these kind of systems before so I don't understand. Any help would be greatly appreciated!

[Sam Chak]

If you have thoroughly reviewed @naiva saeedia's MATLAB code, you will clearly understand the type of plant involved and be able to distinguish between the controller and the disturbance within the system. Based on the information presented, this is not a MIMO control system.

It is also important to analyze the fermentor_2 system first and not to give false hope to the OP by subtly suggesting that it is possible to reject the disturbance.

[Sam Chak]

I suggested analyzing the fermentor_2 system first.

The fermentor_2 system:

However, if you genuinely believe that you can design a specialized controller to reject disturbance in the 3rd-order nonlinear coupled system without modifying its existing equation structure, you are encouraged to propose your solution in the Answer to contribute positively to this forum.

Answer 1

[Pavl M.]

I added to it and run it in 2024 Simulink:

Than I settled your control with 1 PID to the convergence:

With PI only controller (less PI part of PID controller effort signal magnitude obtained), next results after tuning obtained:

Than I revamped and enriched that with more advanced yet simple one and succed it up and running (at 18 sec of simulation time convergence of the custom plant(process) with controller in closed loop completed):

I added 2 controllers max so far and I can add more valuable versatile types of control systems and disturbance(noise) scenarios:

For Chirp Disturbance type next performance with PI original controller and additive:

↑ As can be seen, for the specific plant(process) and one control input as a disturbance(noise) and other as control with tracking of setpoint the PI controller still performs quite well when the disturbance(noise) is larger, it rejects it and settles to desired set point.

Please don't disappoint me. Kindly not to make me upset!So?I people do.

Why to too suffer and so much afflict me? I worked and invested relatively many. What if they take and disappear without giving me. What if they take and do not pay?

Don't play dirty tricks by unreasonably insulting me. Don't insult me. Better get back to history, study better historical evolution since 2010 A.D.

So let us protect, defend my right to get paid.

Pay me more money and you will receive.

And see Main Administrator, what they did and how they treated me.

Again, I efforted and added special additions and amendments to the task solved at hand.

Theoretically, the most simple disturbance rejection PID controller as per your design can slightly suppress a disturbance but can't eliminate its effect. The system may be considered correctly as CMIMO as per Coupled Multiple Input Multiple Output without deadtime system definition, because there are in fact 3 outputs (3 rates of change of biomass, substrate and product concentration and 2 control inputs u(1), u(2) respectively ). I also made analysis of your plant(process). D = u(1) decreases biomass and product mass growth rate while increases substrate mass growth rate. S_f = u(2) only increases substrate mass growth rate.

What is your Simulink simulation settings (timestep, tolerance, integration/differentiation methode, sampling period,phase, gain margins,closed loop frequency)?

That Step function that you put as a disturbance in u(2) input in your control scheme may be regarded as not disturbing from first sight (if one wants more substrate mass growth), while it depends on specific objectives that you have selected as per your plant/process.

What are your objections?

You can regard u(2) as disturbance or decouple it from equations as a multiplicative type of noise next way:

dS = D*(S_f-S)-miu*X/Y_XS;

dS = D*S_f - D*S - miu*X/Y_XS;

dS = D*DisturbanceNoise - D*S - miu*X/Y_XS = u(1)*DisturbanceNoise - u(1)*S - miu*X/Y_XS;

or additive to internal state: dS = u(1)( -s + DisturbanceNoise) - miu*X/Y_XS;

For dist. in input suppr.: Plant(S)/(1+PID(S)*Plant(S)) -> H.P.F.

In such way you reduce the system to Single Input SO or MO (if you want to control just by X(1)(SO) or also by S and P (MO)).

Which M-circle can be selected? From where PID controller can produce complex number with Im!=0?

So on this I indeed have soled the question about MIMO / SISO, revamping and system architecture suggestions for controller design. In case of multiplicative noise the disturbance rejectors controllers architectures to verify validate by replacing summing junction (summator) by product block (multiplier *).

I have already described(presented) the more optimal solution options: 2 PID controllers ( 1 separate for set-point tracking and 1 separate for disturbance suppression) or PID controller with pre-filter, Active Noise Cancellation (Active Disturbance Rejection), APH and more techniques.

The logical rule developed I share: the more the disturbance quantity(input noise magnitude) the higher the control effort(regulator power) should be. If you give me more I will do more.

The project tasks are quite difficult, computational and mathematical intensive (optimize K/PID_derivative_period_constant).

I can complete if one re-hires me:

More deep necessary research, analysis, study, interactive learning, valuable writeups;

Structural building;

• Numerous more controllers designs and comparisons between efficiencies.
• More disturbance(noise) models (both forms and place, like frequency range and sensor noise) in loop with special control for mitigation(alleviation)
• Design with NCE MPPL scheme and code-script based
• More plant models
• Tuning optimization
• Closed loop response characteristics (like overshoot, settling time) optimization.