11. Differential Equations

Specific types of equations with differential terms arise from fundamental relationships such as conservation of mass, energy, and momentum. For example, the accumulation of mass $\frac{dm}{dt}$ in a control volume is equal to the mass inlet $\dot m_{in}$ minus mass outlet $\dot m_{out}$ from that volume.

$\frac{dm}{dt} = \dot m_{in} - \dot m_{out}$

Dynamic models can either be regressed (identified) from data or developed without data with fundamental relationships. Even fundamental relationships may have unknown or uncertain parameters. A combined approach for dynamic system modeling is to use fundamental physics-based relationships with data science. This approach uses the best features of both methods. It creates a model that aligns to measured values and extrapolates to regions where there is no or little data.

The first exercise is to solve differential equations with several examples using odeint. The same examples are also solved with Gekko. For simulation, the two give equivalent results but Gekko is built to use differential equations in optimization or combined with machine learning. The odeint function has a specific purpose to solve ODEs. The function odeint solves the differential equation requires 3 inputs.

y = odeint(model, y0, t)

1. model Function name that returns derivative values at requested y and t values as dydt = model(y,t)
2. y0 Initial conditions of the differential states
3. t Time points where the solution is reported

Solve Differential Equation

We solve the differential equation with initial condition $y(0) = 5$:

$ k \, \frac{dy}{dt} = -y$

where $k=10$. The solution of y is reported from an initial time 0 to final time 20 and with a plot of the result for $y(t)$ versus $t$. Notice how the equation is rearranged to return just the derivative value as dydt = -(1.0/k) * y from the function.

import numpy as np
from scipy.integrate import odeint

# function that returns dy/dt
def model(y,t):
    k = 10.0
    dydt = -(1.0/k) * y
    return dydt

y0 = 5                 # initial condition
t = np.linspace(0,20)  # time points
y = odeint(model,y0,t) # solve ODE

import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(t,y)
plt.xlabel('time'); plt.ylabel('y(t)')
plt.show()

Solve Differential Equation with Gekko

Python Gekko solves the same differential equation problem. Gekko is built for large-scale problems. Additional tutorials on Gekko show how to solve other types of equations and optimization problems.

from gekko import GEKKO

m = GEKKO(remote=False)    # GEKKO model
m.time = np.linspace(0,20) # time points
y = m.Var(5.0); k = 10.0   # GEKKO variable and constant
m.Equation(k*y.dt()+y==0)  # GEKKO equation

m.options.IMODE = 4        # Dynamic simulation
m.solve(disp=False)        # Solve

plt.plot(m.time,y)
plt.xlabel('time'); plt.ylabel('y(t)')
plt.show()

Differential Equation Activity

Solve the differential equation with initial condition $y(0) = 10$:

$ k \, \frac{dy}{dt} = -y$

Compare the five solutions of y from time 0 to 20 with k=[1,2,5,10,20].

Symbolic Solution

Compact differential equation problems may have an analytic solution that can be expressed symbolically. A symbolic math package in Python is sympy. Sympy determines the analytic solution as $y(x)=C_1 \, \exp{\left(-\frac{x}{k}\right)}$. With the initial condition $y(0)=5$, the constant $C_1$ is equal to 5.

from IPython.display import display
import sympy as sym
from sympy.abc import x, k
y = sym.Function('y')
ans = sym.dsolve(sym.Derivative(y(x), x) + y(x)/k, y(x))
display(ans)

Solve Differential Equation with Input u

Differential equations can also have an input (feature) that changes from an external source (exogenous input) such as actively changed by a measurement sensor, a person (manually), or selected by a computer.

Calculate the response y(t) when the input u changes from 0 to 2 at t=5.

$2 \frac{dy(t)}{dt} + y(t) = u(t)$

The initial condition is y(0)=1 and the solution should be calculated until t=15. Tip: The expression y(t) does not mean y multiplied by t. It indicates that y changes with time and is written as a function of time. There are additional examples for ODEINT and GEKKO if you need help.

TCLab Activity

Data Collection

Turn on heater 1 to 100% and record $T_1$ every 5 seconds for 3 minutes. The data should include a total of 37 data points for each temperature sensor.

import numpy as np
import pandas as pd
import tclab
import time
# collect data for 3 minutes, every 5 sec
n = 37
tm = np.linspace(0,180,n)
t1s = np.empty(n); t2s = np.empty(n)
with tclab.TCLab() as lab:
    lab.Q1(100); lab.Q2(0)
    print('Time T1  T2')
    for i in range(n):
        t1s[i] = lab.T1; t2s[i] = lab.T2
        print(tm[i],t1s[i],t2s[i])
        time.sleep(5.0)
# put into dataframe
data = pd.DataFrame(np.column_stack((tm,t1s,t2s)),\
                    columns=['Time','T1','T2'])
data.to_csv('11-data.csv',index=False)

Solve Differential Equation

Use the parameters a, b, and c from 10. Solve Equations or use the following default values.

Parameter	Value
a	78.6
b	-50.3
c	-0.003677

Solve the ordinary differential equation (ODE) with these values.

$\frac{dT_1}{dt} = c (T_1-a)$

The initial condition for $T_1$ is $a + b$. Show the ODE solution at time points between the initial time 0 and the final time 180 sec. Plot the measured $T_1$ on the same plot as the ODE predicted temperature. Add appropriate labels to the plot.