B matrix cannot contain too many zero elements for baltrunc

[JonathanPilgram]

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:

ERROR: LoadError: BoundsError: attempt to access 1-element Vector{BigFloat} at index [1:9]
Stacktrace:
[1] throw_boundserror(A::Vector{BigFloat}, I::Tuple{UnitRange{Int64}})
@ Base ./abstractarray.jl:737
[2] checkbounds
@ ./abstractarray.jl:702 [inlined]
[3] getindex
@ ./array.jl:973 [inlined]
[4] baltrunc(sys::StateSpace{Continuous, BigFloat}; atol::Float64, rtol::Float64, n::Int64, residual::Bool)
@ ControlSystemsBase ~/.julia/packages/ControlSystemsBase/98Jcs/src/matrix_comps.jl:638
[5] top-level scope
@ ~/Documents/Master-AppliedPhysics/MEP/programming/BT.jl:19
in expression starting at /home/jonathan/Documents/Master-AppliedPhysics/MEP/programming/BT.jl:19

The full code is:

using ControlSystems
using HDF5
using GenericLinearAlgebra
setprecision(BigFloat, 256)

data = h5open("test_chain.h5", "r") do file
    A = read(file["A"])
    B = read(file["B"])
    C = read(file["C"])
    D = read(file["D"])
    (A, B, C, D)
end
A = data[1]
A = BigFloat.(A)
B = data[2][1, :]
#B += 0.0000000001 * rand(50, 1) #workaround
C = data[3]
sys = ss(A, B, C, 0)
sys_red, _, _ = baltrunc(sys, n=9) #Line 19
println(sys_red.A)

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?

[baggepinnen]

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

[baggepinnen]

Try calling minreal on the model before baltrunc

[baggepinnen]

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.

[JonathanPilgram]

test_chain.zip
This file contains my data! I will try your suggestions in a bit

[JonathanPilgram]

After calling the sys = sminreal(sys) command before baltrunc I still get the same error...

[baggepinnen]

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]).

[baggepinnen]

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?

[JonathanPilgram]

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.

[baggepinnen]

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
#927

[JonathanPilgram]

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.

[baggepinnen]

I think you can add the branch using

pkg> add ControlSystemsBase#balrealeps

and then restart julia before trying

[baggepinnen]

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

[JonathanPilgram]

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.

[baggepinnen]

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?

[JonathanPilgram]

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:

1. A decay chain model, whereby the initial condition is a certain amount of the first nuclide which then decays along a chain of nuclides. I guess the only real input would be the amount of the first nuclide you choose. However, then I get the errors with which I started, so I chose the input to be an identity matrix.
2. The burnup of a nuclear molten salt reactor. Initially there is a fuel mixture of <10 nuclides. Afterwards the state evolution is made up by decay and by nuclear interactions dependent on neutron flux. For this model I'm not sure how to exactly fit it to control theory. I'm not sure whether to use the initial fuel mixture as input or the neutron flux.

Thanks again for all your help up until now. It means a lot to me :)

[baggepinnen]

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