GeezerVison is showing github.com — Using input_output_response(...), Simple System Gives Erroneous Output · Issue #890 · python-control/python-control
Using input_output_response(...), Simple System Gives Erroneous Output #890DescriptionThis 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()
Activitymurrayrm commented on May 22, 2023 Confirmed this behavior. Using forced_response instead of input_output_response works correctly:
So something odd going on in input_output_response... frohro commented on 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 commented on May 22, 2023 Correct: the input_output_response functions uses scipy.integrate.solve_ivp using its default arguments, which use 'RK45' as the solver and default values for rtol and atol (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. joaoantoniocardoso commented on May 25, 2023 I was asking myself another day -- why not use LSODA as the default? frohro commented on May 25, 2023
MetadataMetadataDevelopmentNo branches or pull requestsIssue actions
|
