TaylorModels and complex evaluation

[lbenet]

This issue is motivated by a discussion which was opened in TaylorSeries, see JuliaDiff/TaylorSeries.jl#288

Quoting @jonniedie:

I have a rational polynomial function G(s) = n(s) / d(s) where n(s) and d(s) are polynomial functions of s. The coefficients of these polynomials have uncertain parameters, which are represented by Intervals. For a given (real) value ω, I would like to evaluate the response G(ω*im) and get the argument and magnitude of the resulting complex number.

When I tried evaluating this directly using Intervals, the resulting bounds were too wide. @mforets recommended I give Taylor models a try, which is when I noticed Complex{<:TaylorModel}s wouldn't work because Complex requires <:Real parameters.

[lbenet]

Here i describe what i think is a work around, using a toy model, which I write as Julia code:

using TaylorModels

a = 1 .. 1;
n(z) = a*z;
d(z) = 1+(a*z)^2;
G(z) = n(z)/d(z);

So we have a parameter a which is an interval, and a function G(z). This function has no singularities in ℝ, so TaylorModel1 should be reliable, though in ℂ things are more subtle. Second, it is easy to compute the real and imaginary parts of G(z), which I write as a two component vector to mimic the real and imaginary parts:

function G(s, ω; a=a)
    aux = 1 + a^2*(s^2-ω^2)
    rr = a*s*(aux + 2*a^2*ω^2)
    ii = a*ω*(aux - 2*a^2*s^2)
    den = 4*a^4*s^2*ω^2 + aux^2
    return [rr/den, ii/den]
end

(I think the algebra and implementation are correct, but do not trust them!)

Now I create TaylorModelNs for the real and imaginary parts:

set_variables("𝓈 ω", order=10);
sm = TaylorModelN(1, 4, IntervalBox(0..0, 2), IntervalBox(-0.5 .. 0.5, 2));
ωm = TaylorModelN(2, 4, IntervalBox(0..0, 2), IntervalBox(-0.5 .. 0.5, 2));

Now, I compute the Taylor model associated to G(s, ω), and then evaluate it on an interval box within the domain:

julia> gm = G(sm, ωm)
2-element Vector{TaylorModelN{2, Interval{Float64}, Float64}}:
  [1, 1] 𝓈 + [-1, -1] 𝓈³ + [3, 3] 𝓈 ω² + [-1.53286, 1.53286]
  [1, 1] ω + [-3, -3] 𝓈² ω + [1, 1] ω³ + [-1.53286, 1.53286]

julia> evaluate(gm, IntervalBox(0..0, 0.125 .. 0.25))
[-1.53286, 1.53286] × [-1.40591, 1.79848]

I think this is the kind of result you are seeking.

[lbenet]

Please let me know if this helps...

[jonniedie]

Thanks @lbenet ! I think I understand most of it, but I may need to play around with it some more to understand how the domain of sm and ωm come into play here. I would expect this example to be exactly zero for the real part of the output, since G in this example can only give purely imaginary outputs for purely imaginary inputs. For that reason, I'm going to alter the toy example slightly to give something that will make real and imaginary parts. I'll also widen the bounds of a a bit, since that's the uncertainty I'm most interested in.

using TaylorModels

a = 0.9 .. 1.1
n(z) = a*z
d(z) = 1 + a*z
G(z) = n(z)/d(z)

To make the evaluation of G a little more general, I'll just quickly make a MyComplex type with just the methods we need to use our existing G

import Base: +, -, *, /, ^

# Not a <:Number for convenience of avoiding ambiguities
struct MyComplex{T}
    re::T
    im::T
end

𝒾 = MyComplex(0, 1)

Base.real(z::MyComplex) = z.re
Base.imag(z::MyComplex) = z.im

Base.to_power_type(z::MyComplex) = z   # For power_by_squaring to work
function Base.inv(z::MyComplex)
    c, d = reim(z)
    den = (c*c + d*d)
    return MyComplex(c/den, -d/den)
end

+(x, z::MyComplex) = MyComplex(x+real(z), imag(z))
*(a, z::MyComplex) = MyComplex(a*real(z), a*imag(z))
*(z::MyComplex, a) = MyComplex(real(z)*a, imag(z)*a)
*(z1::MyComplex, z2::MyComplex) = MyComplex(real(z1)*real(z2) - imag(z1)*imag(z2), real(z1)*imag(z2) + imag(z1)*real(z2))
/(a, z::MyComplex) = a*inv(z)
^(z::MyComplex, n::Integer) = Base.power_by_squaring(z, n)

Base.show(io::IO, z::MyComplex) = print(io, "$(real(z)) + $(imag(z))𝒾")

Now we can make G work as

G(x, y) = [reim(G(x + y*𝒾))...]

The part I'm a little confused about is

set_variables("𝓈 ω", order=10);
sm = TaylorModelN(1, 4, IntervalBox(0..0, 2), IntervalBox(-0.5 .. 0.5, 2));
ωm = TaylorModelN(2, 4, IntervalBox(0..0, 2), IntervalBox(-0.5 .. 0.5, 2));

Do we need such wide bounds on the domains? Or is that just for the example? I was thinking those would just be exact, since the uncertainty we're interested in here is mostly the parametric uncertainty of a. But when I try to set those to some value I want to evaluate at, I seem to be doing something wrong.

julia> p = IntervalBox(0..0, 0.1..0.1);

julia> sm = TaylorModelN(1, 4, p, p)
 [1, 1] 𝓈 + [0, 0]

julia> ωm = TaylorModelN(2, 4, p, p)
 [0.0999999, 0.100001] + [1, 1] ω + [0, 0]

[lbenet]

I'll try to address your comments:

I think I understand most of it, but I may need to play around with it some more to understand how the domain of sm and ωm come into play here.

The variables sm and ωm are the independent TaylorModelNs variables, which represent the real and imaginary parts of the actual complex variable z of the function G(z); this is the core of the work-around for your problem. Each Taylor model consists of a Taylor expansion (in the case of sm and ωm are the independent variables only) around a point (a point in the plane), within a domain (a 2D domain in this case), and a remainder which in this case is an IntervalBox. In the definitions I used for sm and ωm above, the expansion points are given by IntervalBox(0..0, 2), and the domain is IntervalBox(-0.5 .. 0.5, 2). You certainly can adapt these to your problem.

I would expect this example to be exactly zero for the real part of the output, since G in this example can only give purely imaginary outputs for purely imaginary inputs.

You are right, if you have a purely imaginary input G(z) is purely imaginary. But this is not the case I illustrated: The variable sm (real part) is defined around 0, but its domain is -0.5 .. 0.5. So the evaluation of G(sm, ωm) may include non-zero real parts. That's the reason gm[1] is not trivially zero. If you are interested in such a case, you should probably define sm as TaylorModelN(1, 4, IntervalBox(0..0, 0..0), IntervalBox(0..0, -0.5 .. 0.5)); note that the domain in the first component is 0..0.

To make the evaluation of G a little more general,...

I think you the implementation of MyComplex{T} may be useful, but I have to think more about this, though you have to be aware that x + y*𝒾 is the core: if you evaluate sm + im*ωm you will get an error, because the parameter of the Taylor models needs to be real. But your solution may be more general than what I proposed.

Do we need such wide bounds on the domains? Or is that just for the example?

You can change the domain in each direction as well as the expansion point, so the domains may not be identical, nor the expansion points in each direction, and they indeed may be thinner. Note that, while the expansion point may be an interval, it is better it is as thin as possible.

I was thinking those would just be exact, since the uncertainty we're interested in here is mostly the parametric uncertainty of a. But when I try to set those to some value I want to evaluate at, I seem to be doing something wrong.

Since you only care about the parametrical uncertainty of a, you can evaluate G(z) directly, using intervals, which gives you a rigurous result:

julia> G(im*0.1) # using your definition of `G`
[0.00800316, 0.0120028] + [0.088924, 0.109117]im

[lbenet]

I was thinking those would just be exact, since the uncertainty we're interested in here is mostly the parametric uncertainty of a. But when I try to set those to some value I want to evaluate at, I seem to be doing something wrong.

This comment lead me to the question: are you interested in the dependence of a of the results? If so, this may change the way to approach it; I'll try to get back to this later...

[mforets]

@lbenet @jonniedie you work too hard! leave something to discuss in https://juliareach.github.io/juliareach-days-3/ 😄

[jonniedie]

Yes, I'm mostly interested in parametric uncertainty. Doing G(im*0.1) works, but the bounds are too wide, so I was exploring different methods for tightening them. Splicing the interval helped, but I wanted to see what other options there were, because that involved quite a bit of manual work to do (maybe there are more automated ways of doing this that I'm not aware of?).

Also, even though parametric uncertainty is the main driver here, there is some value in starting with an interval domain for ωm because eventually I would like to look at this across a dense range of frequencies split up into sections.

Ultimately the goal is to have some of the functionality (or at least ability to answer the same questions as) the Robust Control Toolbox in MATLAB. ReachabilityAnalysis.jl does a good job of covering robust design and requirements verification for time-domain requirements, but classical linear controls design/analysis is usually done in the frequency domain (or often a mix of the two).

[lbenet]

Thanks @jonniedie for the explanation on the parametric uncertainty.

TL;DR I think TaylorModels can help, but it may get tricky.

I'll use your toy model in the framework I tried to describe initially, which requires using G(s, ω, a) because TaylorModels works with Reals. Let's begin...

using TaylorModels

n(z, a) = a*z
d(z, a) = 1 + a*z
G(z, a) = n(z, a)/d(z, a)

function G(s, ω, a)
    aux = 1 + a*s
    rr = a*(s*aux + a*ω^2)
    ii = a*(ω*aux - a*s*ω)
    den = aux^2+(a*ω)^2
    return [rr/den, ii/den]
end

julia> dom_a = 0.5±0.125 # domain for `a`
[0.375, 0.625]

julia> dom_ω = 1±1/16 # domain for ω
[0.9375, 1.0625]

Note that a is passed as an argument to all functions, because I want to have it as an Interval or TaylorModel-type.

IntervalArithmetic is actually not bad in this case especially (perhaps due to the concrete example), in particular after splicing (mince-ing) the intervals (ω is an interval):

julia> imag( G( dom_ω*im, dom_a) )
[0.243974, 0.591016]

julia> G(0.0, dom_ω, dom_a)[2] # imaginary part, using only reals
[0.243974, 0.591016]

julia> hull(getindex.(G.(0.0, dom_ω, mince(dom_a, 8)),2)) # better
[0.296348, 0.506977]

Now, I'll consider a as a TaylorModel1 variable:

julia> set_taylor1_varname("δa") # use `δa` as the displayed variable

julia> a1 = TaylorModel1(3, 0.5, dom_a )   # independent variable a -> 0.5 + δa
 0.5 + 1.0 δa + [0, 0]

julia> ga1r, ga1i = G(0..0, dom_ω, a1); 

julia> ga1i
 [0.344821, 0.459069] + [0.242038, 0.702564] δa + [-1.16514, -0.273615] δa² + [-0.539255, 1.00459] δa³ + [-0.0376585, 0.037394]

ga1r and ga1i are polynomial approximations of the real and imaginary parts in terms of a, which is expanded around 0.5 and defined in the domain 0.5±0.125. Note that the coefficients are intervals (because
dom_ω is an interval), and they are not specially thin! These expressions can be used as functions, which is nice, in terms of δa:

julia> ga1i(dom_a-mid(dom_a)) # whole domain for `δa` is evaluated
[0.199175, 0.586246]

This result is wider that the one obtained using Intervals. This is because ω is an interval, and therefore all the coefficients of the Taylor expansion become intervals, as I noted above; this is part of the reason that the result with TaylorModels is not so good.

For the specific G(z, a) that we have, let me note that there is only one variable, which is the combination of a*z. So redoing the computation but using that variable (so I'll have to reinterpret the expansion point and the domain):

julia> G1(ra,ia) = [ra, ia] ./ ((1+ra)^2+ia^2);

julia> dom_aω = dom_a * dom_ω # domain of a*ω
[0.351562, 0.664063]

julia> res = G1(0.0, TaylorModel1(3, mid(dom_aω), dom_aω))[2] # mid(dom_aω) is the mid point of dom_aω
 0.40370711824930855 + 0.46903360436370056 δa - 0.699648532292137 δa² + 0.19202799764405293 δa³ + [-0.000369503, 0.000641092]

Note that res is expanded around the mid point of dom_aω, the coefficients of the polynomial are now Float64, and they are inside the interval coefficients of ga1i. In addition, the remainder is much thinner! Then

julia> res(dom_aω-mid(dom_aω))
[0.312237, 0.478368]

This is much better than the naive interval evaluation of G1, and slightly wider than mincing dom_aω in 128 intervals; you can check this. This may be improved by using larger polynomial degree of the TaylorModel1, though I did not check this.

I redo this using both a and ω as TaylorModelNs (we have more than 1 independent variable):

julia> set_variables("δa δω", order=6)
2-element Vector{TaylorN{Float64}}:
  1.0 δa + 𝒪(‖x‖⁷)
  1.0 δω + 𝒪(‖x‖⁷)

julia> aN = TaylorModelN(1, 3, Interval(mid(dom_a)) × Interval(mid(dom_ω)), dom_a×dom_ω)
 [0.5, 0.5] + [1, 1] δa + [0, 0]

julia> ωN = TaylorModelN(2, 3, Interval(mid(dom_a)) × Interval(mid(dom_ω)), dom_a×dom_ω)
 [1, 1] + [1, 1] δω + [0, 0]

julia> gNi = G(0.0, ωN, aN)[2]
 [0.399999, 0.400001] + [0.479999, 0.480001] δa + [0.239999, 0.240001] δω + [-0.704001, -0.703999] δa² + [-0.224001, -0.223999] δa δω + [-0.176001, -0.175999] δω² + [0.179199, 0.179201] δa³ + [-1.13921, -1.13919] δa² δω + [-0.569601, -0.569599] δa δω² + [0.0223999, 0.0224001] δω³ + [-0.000471446, 0.000761742]

julia> gNi( IntervalBox(dom_a-mid(dom_a),dom_ω - mid(dom_ω)) )
[0.309344, 0.479258]

This result is not as good as using the function G1, but better than the naive use of interval arithmetic I described for G above.

[lbenet]

An issue of the approach outlined occurs when G can not be separated in real and imaginary parts easily, which I guess occurs in practice. In that case, your approach with MyComplex{T} seems the way to go.

[jonniedie]

Looks great! Thanks. It looks like the IntervalLinearAlgebra.jl package will help here as well, as sometimes models are given in state-space form, rather than transfer function, and this analysis for those types of models will involve getting the bounding box of the eigenvalues of the state matrix.

So I think this should be enough to get me started on tying these ideas in with ControlSystems.jl to make an interval-based robust control systems package. As soon as I have some free time to work on this, I'll put up a repo and keep everyone posted.