WIP Add transmission-line model

[schillic]

This is based on the paper mentioned in the file and the CORA model and settings.

Problems:

• [x] discretization requires higher Taylor order than in the paper, but CORA implementation also prints Taylor order not high enough when using order 6 (but works with order 8 already)

ERROR: LoadError: AssertionError: c < 1
Stacktrace:
 [1] #_expm_remainder#6(::Int64, ::typeof(IntervalMatrices._expm_remainder),
     ::IntervalMatrix{Float64,IntervalArithmetic.Interval{Float64},
     Array{IntervalArithmetic.Interval{Float64},2}}, ::Float64, ::Int64) at
     .julia/dev/IntervalMatrices/src/exponential.jl:103

• [ ] clarification needed

  • [ ] input factor is ±1 according to paper/CORA implementation; this essentially flips the sign in the plots
  • [ ] X0 is either computed from mid(B) (paper) or from B (CORA implementation)
  • [ ] X0 is either computed with -A⁻¹ (paper) or with +A⁻¹ (CORA implementation)

• [x] initial states are not contained in the flowpipe we found the reason: the implemented algorithm is only sound if the origin is contained in the inputs. so we need the sound version first

  • [x] (likely a consequence:) projection of first set (discretized X0) to U_out (dimension η) already spans interval [-3, 3]*10^4 (paper: [-0.201, 0.201]) the set X0 (called R(0) in the paper) is a zonotope centered in the origin with generator matrix [C 0; 0 D] where C = diag(0.201) and D = diag(1e-3) both in JuliaReach and in CORA so there seems to be a problem in the discretization

  • [x] (likely a consequence:) computed zonotope centers contain NaNs and Infs (most probably just a consequence of the big numbers mentioned above)

• [ ] dimensional errors in projection

further suggestions by @mforets

• change inv(mid(x)) to mid(inv(x))
• change interpretation of X0 (back) to a set interpretation (instead of an interval-matrix interpretation as in CORA)

[mforets]

what does this setting mean?

Seems to me that options.eAt is a cache to hold exp(mid(A) * δ).

computed zonotope centers contain NaNs and Infs

That's bad, that reminds me this issue but i thought it was solved. Do you have an indication where such values arise?

[schillic]

Here is what we want to end up with (CORA results) in the η = 2 (4D) case (corresponding to Figure 9 and 10 of the paper):

For better visibility, here is what happens for time horizon = 0.002 (one time step only):

[schillic]

Here are results with BFFPSV18 and the center matrices:

First time step including the (continuous-time) initial states:

Observations:

• The first plot does not converge to 1 but to -1. This is a consequence of whether we use +1 or -1 as input factor (I checked that) (see the OP).
• Our discretization is much less precise.

[schillic]

Because that precision looked horrible, here is BFFPSV18 with δ' = δ/100.0:

[schillic]

Updated projections of X(0) with improved computation of the exponential remainder (https://github.com/JuliaReach/IntervalMatrices.jl/pull/71):

This is definitely a major improvement but still an order of magnitude away from the true solution.

[schillic]

With the latest fixes, and ignoring a current bug in the discretization for now, the results of ASB07 and ASB07_decomposed (with one-block partition) are equivalent in the first two steps. However, after two steps they diverge, and surprisingly the decomposed version is not a superset of the other.

ASB07

ASB07_decomposed

[schillic]

Another bugfix later, we are getting closer... (blue: ASB07_decomposed, red: ASB07, light green: BFFPSV18 for fixed parameters, yellow: X0)

As can be seen, ASB07_decomposed and ASB07 are very similar until after a while ASB07_decomposed grows (rather quickly). And indeed if we continued the plot, the growth would get arbitrarily big.

Here is the same plot for a 10x smaller time step:

[schillic]

Some updates:

For η = 2 we can make the sets at least not diverge too much by computing the matrix exponential anew in each iteration:

However, η = 4 still diverges:

[schillic]

On the positive side, ASB07_decomposed is much more robust for bigger time steps than BFFPSV18:

δ = 0.0001 (observe that BFFPSV18 and ASB07 produce the same output):

δ = 0.001:

δ = 0.002:

δ = 0.005 (observe the unsoundness between steps 3 and 4 that we still need to fix):

δ = 0.01: