ImperialCollegeLondon/sharpy

scipy used for direct balancing method

AntonioWR opened this issue · 1 comments

Describe the bug
In librom.py under balreal_direct_py, Scipy Schur() decomposition, solve_discrete_lyapunov() and svd() do not accept csr_matrix type for state matrix A.

My temporal way of dealing with it is to transfer it to a dense matrix by the todense() function shown below.

if Schur:
# decompose A
# A is a sparse matrix in csr_matrix format, can not be directly passed into scipy _solver.py so is transformed here first.
A = A.todense()
Atri, U = scalg.schur(A)
# solve Lyapunov
BBtri = np.dot(U.T, np.dot(B, np.dot(B.T, U)))
CCtri = np.dot(U.T, np.dot(C.T, np.dot(C, U)))
Wctri = sollyap(Atri, BBtri)
Wotri = sollyap(Atri.T, CCtri)
# reconstruct Wo,Wc
Wc = np.dot(U, np.dot(Wctri, U.T))
Wo = np.dot(U, np.dot(Wotri, U.T))
else:
# A is a sparse matrix in csr_matrix format, can not be directly passed into scipy _solver.py so is transformed here first.
A = A.todense()
Wc = sollyap(A, np.dot(B, B.T))
Wo = sollyap(A.T, np.dot(C.T, C))

To Reproduce
I am using the Goland wing example, but I change the ROM method to direct balancing method.

If model related, please try and reproduce with one of the provided models (found in /cases/). Else, please provided the input files (a .fem.h5, aero.h5 and case.sharpy).

Expected behavior
A clear and concise description of what you expected to happen.

Screenshots
If applicable, add screenshots to help explain your problem.

System Info (please complete the following information):

  • OS: Ubuntu 20.04
  • SHARPy Version 1.3

Additional context
Add any other context about the problem here.

ngoiz commented

Hi Antonio,

Many thanks for raising this, and apologies for not getting back to you earlier.

I think your approach works well, because I haven't found a scipy alternative to solve the Lyapunov equation for sparse matrices.

If you want you can open a pull request to fix this little issue with your solution. We would greatly appreciate that! (maybe you can add a warning message if the matrices are sparse saying that they need to be dense and that you are transforming them, or something along those lines)

Thanks,
Norberto