This notebook contains material from ND-Pyomo-Cookbook; content is available on Github. The text is released under the CC-BY-NC-ND-4.0 license, and code is released under the MIT license.

# 5.1 Response of a First Order System to Step and Square Wave Inputs¶

Keywords: Simulator, ipopt usage

This notebook demonstrates simulation of a linear first-order system in Pyomo using the Simulator class from Pyomo. The notebook also demonstrates the construction and use of analytical approximations to step functions and square wave inputs.

## 5.1.1 Imports¶

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt

from math import pi

import shutil
import sys
import os.path

if not shutil.which("pyomo"):
!pip install -q pyomo
assert(shutil.which("pyomo"))

if not (shutil.which("ipopt") or os.path.isfile("ipopt")):
!wget -N -q "https://ampl.com/dl/open/ipopt/ipopt-linux64.zip"
!unzip -o -q ipopt-linux64
else:
try:
!conda install -c conda-forge ipopt
except:
pass

assert(shutil.which("ipopt") or os.path.isfile("ipopt"))

from pyomo.environ import *
from pyomo.dae import *

## 5.1.2 First-order differential equation with constant input¶

The following cell simulates the response of a first-order linear model in the form

\begin{align} \tau\frac{dy}{dt} + y & = K u(t) \\ \end{align}

where $\tau$ and $K$ are model parameters, and $u(t)=1$ is an external process input.

In [2]:
tfinal = 10
tau = 1
K = 5

# define u(t)
u = lambda t: 1

# create a model object
model = ConcreteModel()

# define the independent variable
model.t = ContinuousSet(bounds=(0, tfinal))

# define the dependent variables
model.y = Var(model.t)
model.dydt = DerivativeVar(model.y)

# fix the initial value of y
model.y[0].fix(0)

# define the differential equation as a constraint
model.ode = Constraint(model.t, rule=lambda model, t: tau*model.dydt[t] + model.y[t] == K*u(t))

# simulation using scipy integrators
tsim, profiles = Simulator(model, package='scipy').simulate(numpoints=1000)

fig, ax = plt.subplots(1, 1)
ax.plot(tsim, profiles, label='y')
ax.plot(tsim, [u(t) for t in tsim], label='u')
ax.set_xlabel('time / sec')
ax.set_ylabel('response')
ax.set_title('Response of a linear first-order ODE')
ax.legend()
ax.grid(True)

## 5.1.3 Encapsulating into a function¶

In following cells we would like to explore the response of a first order system to changes in parameters and input functions. To facilitate this study, the next cell encapsulates the simulation into a function that can be called with different parameter values and input functions.

In [3]:
def first_order(K=1, tau=1, tfinal=1, u=lambda t: 1):
model = ConcreteModel()
model.t = ContinuousSet(bounds=(0, tfinal))
model.y = Var(model.t)
model.dydt = DerivativeVar(model.y)
model.y[0].fix(0)
model.ode = Constraint(model.t, rule=lambda model, t:
tau*model.dydt[t] + model.y[t] == K*u(t))

tsim, profiles = Simulator(model, package='scipy').simulate(numpoints=1000)

fig, ax = plt.subplots(1, 1)
ax.plot(tsim, profiles, label='y')
ax.plot(tsim, [u(t) for t in tsim], label='u')
ax.set_xlabel('time / sec')
ax.set_ylabel('response')
ax.set_title('Response of a linear first-order ODE')
ax.legend()
ax.grid(True)

first_order(5, 1, 30, sin)

## 5.1.4 Analytical approximation to a step input¶

Keywords: step input

The response to a step change is a common test giving insight into the dynamics of a given system. An infinitely differentiable approximation to a step change is given by the Butterworth function $b_n(t)$

$$b_n(t) = \frac{1}{1 + (\frac{t}{c})^n}$$

where $n$ is the order of a approximation, and $c$ is value of $t$ where the step change occurs.

In [4]:
u = lambda t: 1/(1 + (t/10)**100)

first_order(5, 1, 30, u)
In [5]:
u = lambda t: 1 - 1/(1 + (t/10)**100)

first_order(5, 1, 30, u)

## 5.1.5 Analytical approximation to a square wave input¶

Keywords: square wave input

An analytical approximation to a square wave with frequency $f$ is given by

$$\frac{4}{\pi} \sum_{k=1, 3, 5,\ldots}^N \frac{sin(k\pi/N)}{k\pi/N}\frac{sin(2\pi ft)}{k}$$

where the first term is the Lanczos sigma factor designed to suppress the Gibb's phenomenon associated with Fourier series approximations.

In [6]:
def square(t, f=1, N=31):
return (4/pi)*sum((N*sin(k*pi/N)/k/pi)*sin(2*k*f*pi*t)/k for k in range(1, N+1,2))

u = lambda t: square(t, 0.1)

first_order(5, 1, 30, u)