PRIMA gives an unstable result?

Question

I am working with Modified Nodal Admittance representation of circuits. I am doing Model Order Reduction using PRIMA on MATLAB. I am considering these circuits as Descriptor State-Space systems.

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 :

1. The projection matrix Q : Orthogonal, with the error ~e-15.
2. Checked the poles using MATLAB's pole() : It showed me exactly one unstable pole pair whose Re is ~e+11, so I don't see it crossing over the imaginary axis in any way.
3. 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.
4. hsvd() : shows a couple of unstable modes
5. isstable() : states that its not stable
6. 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 pole() and 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.

Answer 1

I am not able to comment on your answer. However, I ran into similar problem while simulating interconnects in VLSI circuits. That time it turned out I was applying prima incorrectly.

if system is

$E\dot{x}=Ax+Bu$

and Q is projection space, then reduced system is

$Q^{T}EQ\dot{z}=Q^{T}AQz+Q^{T}Bu$

as oppose to

$\dot{z}=Q^{T}E^{-1}AQz+Q^{T}E^{-1}Bu$

Answer 2

The situation is not clear enough. Do you mean the original model is stable? But then it has one unstable pole! If this is the pole of the reduced model, check the eigenvalue of the original model which has smallest real value. It is expected to be close to the imaginary axis. Perhaps the computational error shifts it to the other side of imaginary axis and makes the reduced model unstable. Theoretically, a passive system is a stable system. PRIMA preserves the passivity means it will keep the reduced model stable. It might happen that the original is not (numerically) passive (E=E^T>=0, A^T+A <=0, B^T=C).

Answer 3

After significant muddling through, I figured one important issue.

hsvd() essentially results in the formation of a balanced system. The process of balancing involves taking the chol() of the gramians. chol() requires the input matrix to be positive-definite. Grammians are positive-definite only when the system is controllable(or observable, if you're considering the observability grammian). Being completely controllable and observable means the system is minimal.

My system was most assuredly NOT minimal. So hsvd() was reporting my unobservable/uncontrollable modes as unstable.

I would like others to clarify if this is really a bug. Also refer to my other question, which is not a full follow-up of this question (i.e. I still want an answer to this).