Repository navigation
Using input_output_response(...), Simple System Gives Erroneous Output #890
Description
This system should not have the bump in the output you see about halfway through the simulation.
import matplotlib.pyplot as plt
import numpy as np
import control as ctrl
R = 22e3
C = 470e-6
circuit = ctrl.ss(-1/(R*C), 1/(R*C), 1, 0)
t, resp = ctrl.step_response(circuit)
plt.plot(t, resp)
plt.title('Step Response of the RC Circuit')
plt.grid()
print("The eigenvalue is: ",-1/(R*C))
print("numpy gives the eigenvalue, and eigen vector respectively as: ")
np.linalg.eig(circuit.A)
t = np.arange(0, 120, 0.01)
u = 3.3*148/256-3.3*40/256*np.heaviside(t-60,1)
plt.plot(t,u)
x0 = 3.3*128/256
t, vc = ctrl.input_output_response(circuit, t, u, x0)
plt.plot(t,vc)
plt.grid()
plt.title("Response of RC Circuit")
plt.legend(["u","$v_c$"])
plt.xlabel("time (seconds)")
plt.ylabel("Voltage (Volts)")
plt.show()
Activity
murrayrm commentedon May 22, 2023
Confirmed this behavior. Using
forced_responseinstead ofinput_output_responseworks correctly:So something odd going on in
input_output_response...frohro commentedon May 22, 2023
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.
murrayrm commentedon May 22, 2023
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:
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.
joaoantoniocardoso commentedon May 25, 2023
I was asking myself another day -- why not use LSODA as the default?
frohro commentedon May 25, 2023
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.