Repository navigation
Using input_output_response(...), Simple System Gives Erroneous Output #890
Description
Activity
Confirmed this behavior. Using
forced_responseinstead ofinput_output_responseworks correctly:
So something odd going on ininput_output_response...It appears to me that the input_output_response is using Runga-Kutta 45, and that the tolerance is set to 0.001 or so. I get the same results with that if the rtol is set to low. Here is my python script with the RK45 method. Here is the plot.
`# %% Imports
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp%% Define independent function and derivative function
def x(t):
return 1-np.heaviside(t-60, 1)def f(t, y, c):
R = 22e3
C = 470e-6
dydt = x(t)/R/C - y/R/C
return dydt%% Define time spans, initial values, and constants
tspan = np.linspace(0, 120, 1000)
yinit = [0]
c = []%% Solve differential equation
sol = solve_ivp(lambda t, y: f(t, y, c),
[tspan[0], tspan[-1]], yinit, t_eval=tspan,
rtol = 0.001)%% Plot independent and dependent variable
Note that sol.y[0] is needed to extract a 1-D array
plt.figure(1)
plt.clf()
fig, ax = plt.subplots(num=1)
ax.plot(sol.t, x(sol.t), 'k-', label='Input')
ax.plot(sol.t, sol.y[0], 'k--', label='Output')
ax.legend(loc='best')
plt.show()`With rtol = 1e-5, it looks pretty good.
Correct: the
input_output_responsefunctions usesscipy.integrate.solve_ivpusing its default arguments, which use 'RK45' as the solver and default values forrtolandatol(Scipy default values are 1e-3 for rtol and 1e-6 for atol).The problem you are having is that 'RK45' doesn't detect/handle the discontinuity at t = 60. You can fix this by using 'LSODA' as the method to solve the ODE:
t, vc = ctrl.input_output_response( circuit, t, u, x0, solve_ivp_method='LSODA')which gives the expected response:
It's not clear that the right "fix" is for python-control, but at the least we (I) should add something to the documentation cautioning about the need to be careful with the solver if you have discontinuous inputs.
I was asking myself another day -- why not use LSODA as the default?
It looks like LSODA uses the ratio of the maximum eigenvalue to the minimum eigenvalue to determine the stiffness, and in so this differential equation would not default to a stiff solver with SODA, at least the way I understand it. This problem has difficulty because the forcing function has a sharp discontinuity. Picking a different solver might still solve this though, because the RK45 algorithm is pretty poor for stiff ODE's.
- linked a pull request that will close this issueadd/cleanup documentation on simulation functions #905
on Jun 5, 2023
This system should not have the bump in the output you see about halfway through the simulation.