how can I get the transfer function of this equation using Matlab: dy/dx + K*sqrt(y)-C =0

This is an equation from an inter-acting tanks I used in a closed water circuit microhydro test rig. The first tank is on the ground and has a height of H1. Water is pumped from the first tank. There is a diverting valve that diverts water from the pump back into the first tank at a fraction of the pumped water. The rest of the water from the pump is emtied into the second tank. The height of the water in the second tank is H2. Tank2 it is vertical, and its water exits by gravity into the turbine inside a pipe attached at the bottom of tank2. In detail the system equations are actually:

equation 1: dh1/dt = k1*sqrt(h2) -kp1*(1-x)

equation 2: dh2/dt = kp2*(1-x) -k2*sqrt(1-x)

My problem is to determine if the system is stable or not. In particular, can I use the the transfer function methods to determine the stability of my system?

Answer 1

[Sam Chak]

Hi @Reuel

The square root is a nonlinear term in the system, making the Coupled Tank System nonlinear. Consequently, the transfer function is not defined for a nonlinear system. If you wish to investigate the stability of the system, one approach is to run simulations (Monte Carlo style) over a range of parameters and then graphically determine if the proposed design system is stable.

Additionally, I was unable to interpret your two differential equations with the variable 'x', so I made some modifications.

%% Settings
tspan   = [0 10];     % simulation time
h0      = [1; 0.5];   % initial heights at t = 0

%% Solve IVP
[t, h]  = ode45(@CoupledTank, tspan, h0);

%% Plot results
plot(t, h, 'linewidth', 2), grid on, 
xlabel('t / sec'), ylabel('Water Level / unit')
legend('Water level of Tank 1', 'Water level of Tank 2', 'location', 'e')
title('Water levels in Coupled Tank')

%% Closed-loop Water Circuit Microhydro Test Rig
function dhdt = CoupledTank(t, h)
    % Define State variables
    h1      = h(1);
    h2      = h(2);

    % Parameters
    k1      = 1;
    k2      = 1;
    kp1     = 1;
    kp2     = 1;
    
    % ODE
    dh1dt   = k1*sqrt(h2)  + kp1*(1 - h1);
    dh2dt   = kp2*(1 - h1) + k2*sqrt(1 - h2);
    
    dhdt    = [dh1dt
               dh2dt];
end

[Sam Chak]

Hi @Reuel

If you utilize the linearization method as described by @Anurag, then it can only demonstrate the local stability of the Coupled Tank System around the selected Operating Points.

syms h1 h2 u1 u2 real

% Parameters
k1  = 1;
k2  = 1;
kp1 = 1;
kp2 = 1;

% Coupled-tank system (with dummy external inputs u1 and u2)
f1  = k1*sqrt(h2)  + kp1*(1 - h1)    + u1;
f2  = kp2*(1 - h1) + k2*sqrt(1 - h2) + u2;
F   = [f1
       f2];

% Operating points
h1_op = 1/2*(2 + sqrt(2));
h2_op = 1/2;

% State vector and input vector
X = [h1     % states
     h2];
U = [u1     % inputs
     u2];

% Jacobians for A and B matrices
Al    = jacobian(F, X);
Bl    = jacobian(F, U);
A_lin = subs(Al, {h1, h2}, {h1_op, h2_op});
B_lin = subs(Bl, {h1, h2}, {h1_op, h2_op});

% Linearized State-space model
A   = double(A_lin);
B   = double(B_lin);
C   = eye(2);
D   = 0*C*B;
sys = ss(A, B, C, D)

sys =
 
  A = 
            x1       x2
   x1       -1   0.7071
   x2       -1  -0.7071
 
  B = 
       u1  u2
   x1   1   0
   x2   0   1
 
  C = 
       x1  x2
   y1   1   0
   y2   0   1
 
  D = 
       u1  u2
   y1   0   0
   y2   0   0
 
Continuous-time state-space model.

% Stability analysis
isStable = isstable(sys)    % Can check, but NOT a Proof! Cannot be published in scientific journals

isStable = logical
   1

eigenvalues = eig(A)

eigenvalues = 
  -0.8536 + 0.8280i
  -0.8536 - 0.8280i

pzmap(sys), grid on

Answer 2

[Anurag Ojha]

Hi Reuel

To determine the stability of your system, you can analyze the transfer function of the system. However, before doing that, you need to linearize the nonlinear equations around an operating point. Once you have the linearized equations, you can obtain the transfer function and analyze its stability.

You can refer to following MATLAB code, I have taken certain assumptions while writing the code. Kindly make necessary changes as per your use case.

% Define symbolic variables for states and inputs
syms h1 h2 kp1 kp2 real

% Define hypothetical system equations (replace these with your actual equations)
f1 = -h1 + kp1 * h2; % Example dynamics for h1
f2 = h2 + kp2 * (1 - h1); % Example dynamics for h2

% Define the system equations vector
F = [f1; f2];

% Define the state vector and input vector
X = [h1; h2];
U = [kp1; kp2];

% Calculate the Jacobians for A and B matrices
A = jacobian(F, X);
B = jacobian(F, U);

% Define operating points for linearization 
h1_op = 1;
h2_op = 1;
kp1_op = 0.5;
kp2_op = 0.5;

% Substitute operating points into A and B matrices
A_lin = subs(A, {h1, h2, kp1, kp2}, {h1_op, h2_op, kp1_op, kp2_op});
B_lin = subs(B, {h1, h2, kp1, kp2}, {h1_op, h2_op, kp1_op, kp2_op});

% Convert symbolic matrices to numerical for state-space model
A_num = double(A_lin);
B_num = double(B_lin);

% Assuming direct observation of states and no direct feedthrough
C = eye(2);
D = zeros(2, 2);

% Create the state-space model
sys = ss(A_num, B_num, C, D);

% Display the system for verification
disp('State-space model:');

State-space model:

disp(sys);

  2x2 ss array with properties:

                A: [2x2 double]
                B: [2x2 double]
                C: [2x2 double]
                D: [2x2 double]
                E: []
          Offsets: []
           Scaled: 0
        StateName: {2x1 cell}
        StatePath: {2x1 cell}
        StateUnit: {2x1 cell}
    InternalDelay: [0x1 double]
       InputDelay: [2x1 double]
      OutputDelay: [2x1 double]
        InputName: {2x1 cell}
        InputUnit: {2x1 cell}
       InputGroup: [1x1 struct]
       OutputName: {2x1 cell}
       OutputUnit: {2x1 cell}
      OutputGroup: [1x1 struct]
            Notes: [0x1 string]
         UserData: []
             Name: ''
               Ts: 0
         TimeUnit: 'seconds'
     SamplingGrid: [1x1 struct]

% Analyze the stability and transfer function
isStable = isstable(sys);
tf_sys = tf(sys);
disp(tf_sys);

  2x2 tf array with properties:

       Numerator: {2x2 cell}
     Denominator: {2x2 cell}
        Variable: 's'
         IODelay: [2x2 double]
      InputDelay: [2x1 double]
     OutputDelay: [2x1 double]
       InputName: {2x1 cell}
       InputUnit: {2x1 cell}
      InputGroup: [1x1 struct]
      OutputName: {2x1 cell}
      OutputUnit: {2x1 cell}
     OutputGroup: [1x1 struct]
           Notes: [0x1 string]
        UserData: []
            Name: ''
              Ts: 0
        TimeUnit: 'seconds'
    SamplingGrid: [1x1 struct]

Adding MATLAB documentation for your reference:

• https://www.mathworks.com/help/symbolic/sym.jacobian.html
• https://www.mathworks.com/help/control/ref/tf.html