Closed Jonacool closed 2 months ago
Could you provide the matrix data ro reproduce the error? It will be hard for me to solve it without being able to reproduce it
Try calling minreal
on the model before baltrunc
I think the problem is that you request n=9
modes, but the internal call to balreal
returns only a single mode, i.e., your system may already be non-minimal. minreal
should produce a minimal realization before model-order reduction.
You can also try calling balance_statespace
before calling baltrunc
, which might help with a poorly balanced A
matrix.
test_chain.zip This file contains my data! I will try your suggestions in a bit
After calling the sys = sminreal(sys)
command before baltrunc
I still get the same error...
Your system is structurally non-minimal
julia> sminreal(sys).nx
1
only the first output/state variable is ever affected by the input. Why do you set B
to the initial condition of a simulation? Are you trying to simulate an impulse response? If so, set B = I
, perform the balanced truncation and then call impulse(sys_red[:, 1])
.
I'll close this issue since I don't believe there's anything wrong in the implementation. #926 adds a more friendly error message.
Please feel free to continue the discussion here until you have found a solution to your original problem. To help you, a bit more context would be useful, what is your final goal, what does the model represent and where is it coming from?
Thanks for you help!
I'm trying to make a reduced order model of a fuel burnup of a nuclear reactor. Before I tried the proper orthogonal decomposition, however the timescales in a nuclear reactor vary wildly (10^-20 - 10^40), so it regularly creates unstable models. That is why I'm trying balanced truncation in Julia, with BigFloats. My goal with the decay chain is just to test the methodology, so I'm mainly interested in the impulse response, however in general I would need to make a model that can handle all kinds of initial conditions.
If I set B to I, the accuracy of the model seems to suffer and it still gives a similar error to the original error if I set the difference between the decay constants large enough (40 orders of magnitude), only it says that it is trying to access an 8 element vector.
Do you need to reduce the complete model all at once, or could you keep your time scales separate and reduce the fast and slow parts of the model separately before combining?
You can try this PR branch and see if you can retain more modes than 8 https://github.com/JuliaControl/ControlSystems.jl/pull/927
How do I use the PR branch?
But no, I can't separate the time scales, sometimes a slow decaying nuclide comes from a fast decay and vise-versa.
I think you can add the branch using
pkg> add ControlSystemsBase#balrealeps
and then restart julia before trying
Another thing you can try is
using RobustAndOptimalControl
sys_fast, sys_slow = stab_unstab(sys, smarg = -100);
this splits the system additively into two systems with different time scales, one with poles faster than 100 (you may need to tweak this value) and one with slower. You could then reduce the order of the two models independently, and later combine them again as
sys_red = sys_fast_red + sys_slow_red
Something like this
s1,s2 = stab_unstab(sys, smarg = -1e4);
# Figure out a reasonable time scale to split the systems at
s1r, _ = balreal(s1);
s2r, _ = balreal(s2);
@show s1.nx, s2.nx
@show s1r.nx, s2r.nx
# smarg = -1e4 seems reasonable, now reduce fast and slow dynamics independently
s1r, _ = baltrunc(s1, n=5);
s2r, _ = baltrunc(s2, n=5);
sys_red = s1r + s2r # Final reduced-order system of order 5+5
Thanks for the suggestions again. I finally got around to implementing them.
Splitting at -1e15
gave me errors that the balanced realization has not enough dimensions left after making a balanced realization for the fast part. I set the eigenvalues of the original matrix to be logarithmically spread between [-1;-1e30]
---I chose the splitting to be halfway between these values, as that seems sensible. Maybe it has something to do with how the splitting is calculated and the difference between the eigenvalues (stiffness) being so large?
Your other solution gives a ROM where the initial condition decays after a single time and flatlines.
I think the crux of the problem is that the time scales vary much more than the machine precision, so many of the calculations are not precise and sometimes the system is unstable, when expressed in machine precision, whilst system should not be unstable.
Did you perform the computations above using BigFloats?
The model-reduction methods in this package all assume an input-output model, but this does not appear to fit your use case, you don't have any inputs?
The simulation is done in Python float64 using the scipy ivp solver with BDF method. All the calculations in Julia are done using BigFloats.
I'm not entirely sure what my inputs must be, as I'm quite new to control theory. I have two models I'm trying to simulate:
Thanks again for all your help up until now. It means a lot to me :)
The simulation is done in Python float64
You could perform the simulation in BigFloat in julia to see if that resolves some of these issues.
The models you describe are not input-output models, they are autonomous ODEs (no input). The methods in this package mostly assume IO models, in particular, the methods for model reduciton. You can simulate autonomous models just fine though, in which case the initial condition is provided as initial condition and not as input. For such a model, the input matrix is empty zeros(nx, 0)
and the initial condition is provided through the keyword argument x0
to the function lsim
During balanced truncation of a decay chain model I found certain errors depending on my choice of B. If it contains to many zeros, then I get an bounds error:
The full code is:
As I generate A,B,C using python I will describe them. A is a 50x50 sparse matrix with negatives on the diagonal and same values, but positives just right next to it. B should be the initial condition, which is only element 1 being 1 and the other 49 zero. C is a 50x50 identity matrix.
I found a workaround, by adding a tiny value to every element of B, but I do not feel comfortable doing that. Why does the code 'need' the nonzero elements?