Discussion of frequency responses

[olof3]

@baggepinnen @mfalt and I have been starting to look at how to utilize the specialized inversion algorithm for Hessenberg matrices that was in introduced in 1.3 to speed up the computation of frequency responses. Here are some thoughts to start of further discussions.

A question is how to exploit properly and avoid re-compuing Hessenberg factorziations too often.

1. freqresp should obvioulsy exploit the Hessenberg structure if the StateSpace is known to be structured (i.e., specific subtypes of HeteroStateSpace). However what should be done if A is a regular matrix?
   A. Should the user have to explicitly have construct a structured statespace with a Hessenberg A matrix?
   B. Should the Hessenberg factorization somehow be cached in the StateSpace object?
   C. Automatically compute a Hessenberg factorization for each call to freqresp and evalfr (probably overloading the broadcase evalfr. for StateSpace). Either a strucatured statespace could be constructed, or the factorization could just be used directly.
   D. Other ideas?

To me, C seems like the best option, it is much more straightforward than B, and if the user really cares about efficiency it should be relatively easy to compute a structured statespace with Hessenberg structure.

Some other thoughts while we are discussing frequency responses

2. What ordering of the output of freqresp dimensions should we use? To me it seems seems more natural to have the index corresponding to frequency first, as typically one considers the frequencies one at a time. (except when plotting, which is why we decided on this ordering in the first place, I guess)

3. It should be possible to use freqresp on a single number (with a matrix output). This situation occurs quite often to me and it is a hassle having to go through evalfr or dropping dimensions.

4. The name of evalfr is notoriously easy to confuse with freqresp. Some ideas: i) rename it to evaltf, or something even better (evaltf has some obvious problems, but I'd say it is much better than evalfr). ii) Only allow Complex inputs in the exported function (we already have dcgain for the by far most common real input. Or is there something like DualNumbers for which this would break?

[baggepinnen]

To me, C seems like the best option, it is much more straightforward than B, and if the user really cares about efficiency it should be relatively easy to compute a structured statespace with Hessenberg structure.

1. Right now the approach is to automatically compute (but not yet make use of) the Hessenberg factorization when calling freqresp and noting in the docstring for evalfr that it's more efficient to use freqresp for vectors of frequencies. I think this approach is okay if we make good use of the Hessenberg factorization in freqresp.

2. As long as it's the same conventions as is used by all other function, like bode, nyquist etc. The distinction only matters for !SISO systems, for SISO one can just call vec (or make use of autovec [WIP] Autovec #209 ). What does one typically want to do with the result of freqresp?

• Plot it
• Reduce it to a scalar, in cost function or similar
  The question is really in situations that are performance critical, do you operate per frequency or per "subsystem" (input-output pair)?

3. This is kind-of what evalfr does, but evalfr require a complex frequency.

4. I agree that evalfr is a footgun. I do not think it's a good idea to restrict to imaginary inputs in evalfr since I might be in a position where I want to evaluate it for a real number, but changing the name is a good idea. I find evaltf an okay name in lack of a better one. Dual is notoriously hard to combine with Complex, but I'm not sure if we get around that either way? Maybe with the proper use of cis(x) instead of exp(im*x) everywhere?

[olof3]

1. Right now the approach is to automatically compute (but not yet make use of) the Hessenberg factorization when calling freqresp and noting in the docstring for evalfr that it's more efficient to use freqresp for vectors of frequencies. I think this approach is okay if we make good use of the Hessenberg factorization in freqresp.

Good.

2. As long as it's the same conventions as is used by all other function, like bode, nyquist etc. The distinction only matters for !SISO systems, for SISO one can just call vec (or make use of autovec [WIP] Autovec #209 ). What does one typically want to do with the result of freqresp?

   • Plot it

   • Reduce it to a scalar, in cost function or similar
     The question is really in situations that are performance critical, do you operate per frequency or per "subsystem" (input-output pair)?

Adding to the list

• Computing the singular values of the transfer function matrix at every frequency
• Frequency-wise multiplyication with other frequency responses
• Analyzing the frequency response of individual channels

I think that the frequency-by-frequency format is to be preferred for performance critical operations, but I guess this is only one aspect. The print way it prints and what is most natural are also quite important.

It would be nice if frequency responses and time responses were handled consistently.

3. This is kind-of what evalfr does, but evalfr require a complex frequency.

Having to multiply with a complex number is somewhat annoying but for discrete-time systems it is really quite a headache, then it is almost easier with freqresp(G, [w])[1]. I also don't see any good reason for not overloading freqresp.

4. I agree that evalfr is a footgun. I do not think it's a good idea to restrict to imaginary inputs in evalfr since I might be in a position where I want to evaluate it for a real number, but changing the name is a good idea. I find evaltf an okay name in lack of a better one. Dual is notoriously hard to combine with Complex, but I'm not sure if we get around that either way? Maybe with the proper use of cis(x) instead of exp(im*x) everywhere?

Good, yes, evaltf would be decent, but we can think some more about it. If we do a name change there would be little need to worry about restricting to complex etc.

[baggepinnen]

Crossref #261
@olof3 maybe the contents of that issue can be reproduced here and the old issue closed?

[olof3]

Yes, this issue seems to deal with everything except

Remove the keyword argument in (sys::TransferFunction)(z_or_omega::Number, map_to_unit_circle::Bool), the same thing can be achieved by the proposed method for freqresp.

So let's close #261

[baggepinnen]

Note to self, time to make use of that hessenberg stuff..

Our freqresp is extremely slow for large systems

julia> @time mag1,_ = bode(G, w);
  8.795737 seconds (97.19 k allocations: 16.532 GiB, 12.24% gc time)

julia> @time magd = abs.(DescriptorSystems.freqresp(Gd, w));
  0.084311 seconds (6.31 k allocations: 7.584 MiB)

julia> G.ny, G.nx, G.nu
(3, 270, 3)

[baggepinnen]

freqresp now uses hessenberg

• (actually) use hessenberg factorization for statespace freqresp #598

the output convention is the same as before, but the output comes in an PermutedDimsArray which is a non-allocating permutation of the layout ny, nu, nw. An ny, nu, nw array can be ontained by accessing the .parent field.

One can also call freqresp on a single number and get a matrix.