scipy / scipy

SciPy library main repository
https://scipy.org
BSD 3-Clause "New" or "Revised" License
13.16k stars 5.21k forks source link

BUG: dstep doesn't work with dlti with complex coefficients. #18273

Open fernandohds564 opened 1 year ago

fernandohds564 commented 1 year ago

Describe your issue.

When a dlti system with complex coefficients is created, the functions dstep and the method step from the dlti object don't work properly. I think it is because of a bug in function dlsim, which assumes real arrays for the object states.

The code below creates this image, where we can see that when we have a complex an discrete system the imaginary part of the step-response is zero and the real part does not match the ones generated with equivalent discrete real transfer function and the continuous transfer function with complex coefficients.

image

Reproducing Code Example

import scipy.signal as scysig
import matplotlib.pyplot as mplt
import numpy as np

w0 = 2
damp = 0.5 * w0
wr = np.sqrt(w0**2 - damp**2)

G = scysig.lti([damp], [1, +damp - 1j*wr])
Gr = scysig.lti([damp, damp**2], [1, 2*damp, w0**2])
Gi = scysig.lti([damp*wr], [1, 2*damp, w0**2])
Gd = G.to_discrete(dt=0.1, method='zoh')
Grd = Gr.to_discrete(dt=0.1, method='zoh') 
Gid = Gi.to_discrete(dt=0.1, method='zoh') 

t, r = G.step()
_, rr = Gr.step(T=t)
_, ri = Gi.step(T=t)
# td, (rd, ) = scysig.dstep(Gd, t=t)
td, (rd, ) = Gd.step(t=t)
_, (rrd, ) = Grd.step(t=t)
_, (rid, ) = Gid.step(t=t)
fig, (ax, ay) = mplt.subplots(2, 1, sharex=True)

ax.plot(t, r.real, label='Continuous Complex')
ax.plot(t, rr.real, '--', label='Continuous Real')
ax.plot(td, rd.real, label='Discrete Complex')
ax.plot(td, rrd.real, '.', label='Discrete Real')
ax.set_ylabel('Real Response')
ax.grid(True, alpha=0.5, ls='--')
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize='small')

ay.plot(t, r.imag, label='Continuous Complex')
ay.plot(t, ri.real, '--', label='Continuous Imag.')
ay.plot(td, rd.imag, label='Discrete Complex')
ay.plot(td, rid.real, '.', label='Discrete Imag.')
ay.set_ylabel('Imaginary Response')
ay.set_xlabel('Time [s]')
ay.grid(True, alpha=0.5, ls='--')
ay.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize='small')

fig.tight_layout()
fig.show()

Error message

/opt/mamba_files/mamba/envs/sirius/lib/python3.9/site-packages/scipy/signal/_filter_design.py:1746: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless
  warnings.warn("Badly conditioned filter coefficients (numerator): the "
/opt/mamba_files/mamba/envs/sirius/lib/python3.9/site-packages/scipy/signal/_ltisys.py:3512: ComplexWarning: Casting complex values to real discards the imaginary part
  xout[i+1, :] = (np.dot(system.A, xout[i, :]) +
/opt/mamba_files/mamba/envs/sirius/lib/python3.9/site-packages/scipy/signal/_ltisys.py:3514: ComplexWarning: Casting complex values to real discards the imaginary part
  yout[i, :] = (np.dot(system.C, xout[i, :]) +
/opt/mamba_files/mamba/envs/sirius/lib/python3.9/site-packages/scipy/signal/_ltisys.py:3518: ComplexWarning: Casting complex values to real discards the imaginary part
  yout[out_samples-1, :] = (np.dot(system.C, xout[out_samples-1, :]) +

SciPy/NumPy/Python version and system information

1.10.0 1.23.5 sys.version_info(major=3, minor=9, micro=2, releaselevel='final', serial=0)
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
  NOT AVAILABLE
openblas_clapack_info:
  NOT AVAILABLE
flame_info:
  NOT AVAILABLE
atlas_3_10_threads_info:
  NOT AVAILABLE
atlas_3_10_info:
  NOT AVAILABLE
atlas_threads_info:
  NOT AVAILABLE
atlas_info:
  NOT AVAILABLE
lapack_info:
    libraries = ['lapack', 'blas', 'lapack', 'blas']
    library_dirs = ['/opt/mamba_files/mamba/envs/sirius/lib']
    language = f77
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
  NOT AVAILABLE
atlas_3_10_blas_threads_info:
  NOT AVAILABLE
atlas_3_10_blas_info:
  NOT AVAILABLE
atlas_blas_threads_info:
  NOT AVAILABLE
atlas_blas_info:
  NOT AVAILABLE
blas_info:
    libraries = ['cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['/opt/mamba_files/mamba/envs/sirius/lib']
    include_dirs = ['/opt/mamba_files/mamba/envs/sirius/include']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    libraries = ['cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['/opt/mamba_files/mamba/envs/sirius/lib']
    include_dirs = ['/opt/mamba_files/mamba/envs/sirius/include']
    language = c
lapack_opt_info:
    libraries = ['lapack', 'blas', 'lapack', 'blas', 'cblas', 'blas', 'cblas', 'blas']
    library_dirs = ['/opt/mamba_files/mamba/envs/sirius/lib']
    language = c
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    include_dirs = ['/opt/mamba_files/mamba/envs/sirius/include']
ilayn commented 1 year ago

Linear systems with complex coefficients are not supported.

tupui commented 1 year ago

@ilayn is it worth adding a note about complex coefficients in the documentation? Also is this something we want to keep as an enhancement?

fernandohds564 commented 1 year ago

Hi guys. Sorry, I didn't know that complex coefficients were not supported. Would you prefer that I close this issue?

Just as a feedback from a user point of view, I think it would be useful to have support for complex coefficients in the future. While I don't think they're strictly necessary, being able to use them many times makes modeling simpler and more intuitive.

Thank you for the quick reply, and all the work developing the library!