by Jeffrey Kantor (jeff at nd.edu). The latest version of this notebook is available at https://github.com/jckantor/CBE30338.
In the example we show how to fit the step response of a nonlinear system, a gravity drained tank, to a first order linear system.
For a tank with constant cross-sectional area, such as a cylindrical or rectangular tank, the liquid height is described by a differential equation
$$A\frac{dh}{dt} = q_{in}(t) - q_{out}(t)$$where $q_{out}$ is a function of liquid height. Torricelli's law tells the outlet flow from the tank is proportional to square root of the liquid height
$$ q_{out}(h) = C_v\sqrt{h} $$Dividing by area we obtain a nonlinear ordinary differential equation
$$ \frac{dh}{dt} = - \frac{C_V}{A}\sqrt{h} + \frac{1}{A}q_{in}(t) $$in our standard form where the LHS derivative appears with a constant coefficient of 1.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
Cv = 0.1 # Outlet valve constant [cubic meters/min/meter^1/2]
A = 1.0 # Tank area [meter^2]
# inlet flow rate in cubic meters/min
def qin(t):
return 0.15
def deriv(h,t):
return qin(t)/A - Cv*np.sqrt(h)/A
IC = [0.0]
t = np.linspace(0,200,101)
h = odeint(deriv,IC,t)
plt.plot(t,h)
plt.xlabel('Time [min]')
plt.ylabel('Level [m]')
plt.grid();
The step response of the gravity drained to a change in flowrate looks similar to the step response of a firat order linear system. Let's try a linear approximation
$$\tau\frac{dx}{dt} + x = Ku$$which has a step response solution that can be written
$$x(t) = x_{ss} + (x_0 - x_{ss})\exp(-t/\tau)$$where $x_{ss} = Ku_{ss}$. There are two parameters, $K$ and $\tau$, which we need to estimate in order to fit the linear approximation to the nonlinear simulation results computed above.
The steady state gain $K$ of the linear system is given by
$$ K = \frac{x_{ss} - x(0)}{u_{ss} - u_0}$$where $u_0$ is the initial input, $u_{ss}$ is the steady-state input, and $x_0$ and $x_{ss}$ are corresponding values of the state variable. In the case of liquid level, $h\sim x$ and $q_{in}\sim u$, therefore an estimate of $K$ is
$$ K = \frac{h_{ss} - h_0}{q_{in,ss} - q_{in,0}}$$q0,h0 = 0,0 # initial conditions
qss = qin(t[-1]) # final input
hss = h[-1] # python way to get the last element in a list
K = (hss-h0)/(qss-q0) # step change in output divided by step change in input
print('Steady-State Gain is approximately = ', K)
Steady-State Gain is approximately = [ 14.98595148]
From the formula for the solution of a first-order linear equation with constant input,
$$\frac{x_{ss} - x(t)}{x_{ss} - x_0} = \exp(-t/\tau) \qquad \implies \qquad \tau = \frac{-t}{\ln\frac{x_{ss} - x(t)}{x_{ss} - x_0}}$$We pick one point representative of the transient portion of the nonlinear response. In this case the response at $t = 25$ minutes accounts for $\approx$60% of the ultimate response, so we choose point as a representative point.
k = sum(t<25) # find index in t corresponding to 25 minutes
tk = t[k]
hk = h[k]
tau = -tk/np.log((hss-hk)/(hss-h0))
print('Estimated time constant is ', tau)
Estimated time constant is [ 23.92271579]
u0 = q0
uss = qss
xss = K*(uss - u0)
xpred = xss - xss*np.exp(-t/tau)
plt.plot(t,h)
plt.plot(t,xpred)
plt.legend(['Nonlinear Simulation','Linear Approximation'])
plt.xlabel('Time [min]')
plt.ylabel('Level [m]')
plt.title('Nonlinear Simulation vs. Linear Approximation');
