Azzera filtri
Azzera filtri

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

18 visualizzazioni (ultimi 30 giorni)
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?

Risposta accettata

Sam Chak
Sam Chak il 12 Giu 2024
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
  1 Commento
Sam Chak
Sam Chak il 12 Giu 2024
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

Accedi per commentare.

Più risposte (1)

Anurag Ojha
Anurag Ojha il 12 Giu 2024
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:

Prodotti


Release

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by