bug: delayed dirac input not considered in forced_response

Hi,

I have encountered a strange issue using the function forced_response with a NonlinearIOSystem object:

When applying a dirac as input, if the dirac is applied after too much time after the beginning of the simulation, the response of the system does not include its impact.
Here is a small piece of code to reproduce the issue:

import numpy as np
import matplotlib.pyplot as plt
import control

def dynamic_function(t, x, u, param=None):
    """
    Dynamic function for the system.
    """
    A = np.array([[0, 1], [-1, -1]])
    B = np.array([[0], [1]])
    return A @ x + B @ u

def output_function(t, x, u, param=None):
    """
    Output function for the system.
    """
    C = np.array([[1, 0]])
    return C @ x

sys = control.NonlinearIOSystem(
    dynamic_function, output_function,
    inputs=1,
    outputs=1,
    states=2
    )

t = np.linspace(0, 10, 100)
u = np.zeros_like(t)


dirac_index = 6


u[dirac_index] = 1  # Dirac input at t=0

results = control.forced_response(sys, T=t, U=u)

plt.figure()
plt.plot(t, results.y[0], label='Output')
plt.title('Response of Nonlinear System to Dirac Input')
plt.xlabel('Time (s)')
plt.ylabel('Output')
plt.grid()
plt.legend()
plt.show()

If dirac_index is inferior or equal to 5, the normal output is obtained while if dirac_index is superior to 6, both the output and the state of the function are null for the whole simulation time.

Here are the installed package:
├── control v0.10.1
│ ├── matplotlib v3.10.3
│ ├── numpy v2.2.5
│ └── scipy v1.15.3

Using python 3.12.3.

I took a look at this. The forced_response function calls input_output_response for a nonlinear system, which in turn calls scipy.integrate.solve_ivp. It is that function that is "ignoring" the input. As noted in the python-control documentation:

If discontinuous inputs are given, the underlying SciPy numerical integration algorithms can sometimes produce erroneous results due to the default tolerances that are used. The ivp_method and ivp_keywords parameters can be used to tune the ODE solver and produce better results. In particular, using 'LSODA' as the ivp_method or setting the rtol parameter to a smaller value (e.g. using ivp_kwargs={'rtol': 1e-4}) can provide more accurate results.

Re-running the example and following these directions fixes this specific issue for the case of dirac_index = 6:

u = np.zeros_like(t)
dirac_index = 6
u[dirac_index] = 1

resp = control.input_output_response(
    sys, T=t, U=u, solve_ivp_method='LSODA')
resp.plot(label='delay=6')

However, the simulation still fails if you use a longer delay (I tried 27) and, interestingly, also fails if you use a delay of 5:

Based on this, it appears the problem is not in python-control per se, but rather in solve_ivp. I suspect that the issue is that the solve_ivp methods are using adaptive step sizes and since the system doesn't change initially, the step size gets large and "skips over" the input. To test that, I used the max_step keyword to limit the maximum step size to 0.1 (the width of the input pulse):

u = np.zeros_like(t)
dirac_index = 27
u[dirac_index] = 1

resp = control.input_output_response(
    sys, T=t, U=u, solve_ivp_kwargs={'max_step': 0.1})
resp.plot(label='delay=27, max_step=0.1')

This generates the correct plot:

I'll create a small PR to update the documentation to add this information (and fix up some typos that I found while digging into this).

Thanks a lot for this complete answer!

update documentation for input_output_response to address gh1152 (pyt…

Navigation Menu

Repository navigation

bug: delayed dirac input not considered in forced_response #1152

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

bug: delayed dirac input not considered in forced_response #1152

Description

Activity

murrayrm commented on Jun 22, 2025

BobAubouin commented on Jun 25, 2025

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions