I have been working with this for quite a while. It is a guaranteed passivity preserving algorithm.
However, for some circuit systems, the implementation I have misbehaves. Although I pass in a stable model, it returns an unstable model(!). So I tested the input and output in multiple ways.
Things I have checked :
- The projection matrix Q : Orthogonal, with the error ~e-15.
- Checked the poles using MATLAB's
pole(): It showed me exactly one unstable pole pair whose
Reis ~e+11, so I don't see it crossing over the imaginary axis in any way.
- Checked the poles using
eig(-G\C)— the poles are actually the inverse of these values, i.e. of those which are non-zero : shows a SINGLE pole (not pair) at ~e-19.
hsvd(): shows a couple of unstable modes
isstable(): states that its not stable
- A custom checking function (from my senior student guide) : 1 unstable pole.
The poles of the original system are also ~e+11. The original is absolutely stable (and passive I think (transmission line), but that's just my speculation.)
Now I am really confused. What is the difference between
eig(-G\C)? And that one pole according to the latter seems quite clearly the real problem.
E(or C) is symmetric positive definite, and A(or G) satisfies G+G^T to ~e-15. However, B != C^T, the algorithm I'm working on is specifically targeted towards inputs != outputs scenario. I have repeatedly enquired and searched, but haven't seen anything about this condition being necessary – the authors just assumed simply this condition for simplicity, as far as I know.
I hope someone can give me further suggestions on what to check for.
Edit: I find that this is happening for only some values of reduced order model. For others,
isstable() reports unstable and
eig(-G\C) says all clear. For yet others,
isstable() reports stable, but
hsvd() still finds unstable modes.
I don't expect people to help me debug this exact issue (that would be asking too much), but I want some help with identifying known gotchas, and discovering unknown ones.