System idenfitifaction of data in the frequency domain

Question

System idenfitifaction of data in the frequency domain

KronosTheLate opened this issue 4 years ago · 14 comments

I have a dataset, which contains the data of a bodeplot of some system. So I have a bunch of frequencies, along with the corresponding amplitudes and phases. I want to use this package to estimate it's transfer function. I do however only see functions in this package, and in the examples, for estimating transfer functions (or StateSpace models, which should be same-same as they are easily converted) for time-signals. Is there functionality for estimation based on frequency responses? Or do anyone have suggestions for other ways to estimate a system for my data?

I can think of the possibility of doing an inverse Fourier transform to get my data over to a time-signal, but that seems like more work than what should be necessary. Maybe, if this functionality is not implemented, this would be a way to add the it?

Answer 1 · 2020-10-23T16:22:31.000Z

A bode curve is a transfer function. so I'm going to assume that you are looking for a parametric representation of the system, such as a low-order rational transfer function.

I currently have no method in this package for the identification you describe. In general, this is a non-convex problem that is typically solved using a gradient-based method like LM. I do some identification like this in a different package (not public) and it can indeed be quite useful, so maybe it would be a good idea to at least provide som rudimentary functionality here.

I have never tried to inverse transform the data to time domain and perform the identification there, I guess it would work though, at least to form the initial guess for a gradient-based refinement.

Answer 2 · 2020-10-23T19:21:26.000Z

A bode curve is a transfer function. so I'm going to assume that you are looking for a parametric representation of the system, such as a low-order rational transfer function.

I am indeed looking for a low order rational transfer function (1 zero 2 poles).

I currently have no method in this package for the identification you describe. In general, this is a non-convex problem that is typically solved using a gradient-based method like LM. I do some identification like this in a different package (not public) and it can indeed be quite useful, so maybe it would be a good idea to at least provide som rudimentary functionality here.

I don't know what LM is, but I would love some rudimentary functionality for finding a transfer function from my signal. It is my assignedment in a class, and as MatLab can do it, I want to be able to "keep up" in Julia, and show/know that Julia can do all the things asked of us that MatLab can.

I have never tried to inverse transform the data to time domain and perform the identification there, I guess it would work though, at least to form the initial guess for a gradient-based refinement.
I realize that calling it Fourier was completely wrong. I just had a 4 hour lecture on it just before writing this, whoops.
But yhea, making the transformation to the time-domain does seem like a possibility, but not the most elegant of solutions.

So do you think that you could, without too much hassle, add your private package functionality for doing what I am looking for? Because that sounds great for me at least ^_^

Answer 3 · 2020-10-23T20:04:31.000Z

It'll take me some time, the private code is not suitable for sharing. It's not too much to setup though, the following (untested) should be close to what's needed

using Optim, ComponentArrays, ControlSystems, ControlSystemIndentification

function loss(p::ComponentVector, data::FRD)
    a,b = p.a, p.b
    G = tf(b,a)
    mag = vec(freqresp(G, data.w))
    @. mag = log(abs(mag)) - log(abs(data.r))
    return sum(abs2, mag)
end

p0 = ComponentVector((a=[1.0,1.0,1.0], b = [1.0]))

opts = Optim.Options(store_trace=true, show_trace=true, show_every=1, iterations=10, allow_f_increases=false, time_limit=100, x_tol=0, f_tol=0, g_tol=1e-8, f_calls_limit=0, g_calls_limit=0)

Optim.optimize(p->loss(p,data), p0, BFGS(), opts, autodiff=:forward)

Answer 4 · 2020-10-23T20:09:30.000Z

Yeah, it works well

using Optim, ComponentArrays, ControlSystems, ControlSystemIdentification

function loss(p::ComponentVector, data::FRD)
    a,b = p.a, p.b
    G = tf(b,a)
    mag = vec(freqresp(G,data.w))
    @. mag = log(abs(mag)) - log(abs(data.r))
    return sum(abs2, mag)
end

# Generate fake data
w = exp10.(LinRange(-1, 2, 100))
Gtest = tf([1.2], [1.1, 0.9, 1.2])
data = FRD(w, Gtest)

p0 = ComponentVector((a=[1.0,1.0,1.0], b = [1.0])) # Initial guess

opts = Optim.Options(store_trace=true, show_trace=true, show_every=1, iterations=10, allow_f_increases=false, time_limit=100, x_tol=0, f_tol=0, g_tol=1e-8, f_calls_limit=0, g_calls_limit=0)

res = Optim.optimize(p->loss(p,data), p0, BFGS(), opts, autodiff=:forward)

Gest = tf(res.minimizer.b, res.minimizer.a)


pole(Gtest)
pole(Gest)

bodeplot(Gtest, w)
bodeplot!(Gest, w)

Answer 5 · 2020-10-24T05:20:25.000Z

It's implemented now, see
https://github.com/baggepinnen/ControlSystemIdentification.jl#parametric-estimation
Feel free to try it out and let me know if it works.

Answer 6 · 2020-10-24T09:11:54.000Z

Feel free to try it out and let me know if it works.

I had a few problems with precompilation, one of which was the repeated printing of the error seen in the picture below:

But after shutting down the process and trying multiple times, it precompiled.

Other than that problem, there are some things I don't understand in the notebook example. I have to provide a sample time, but as my data is in the frequency domain, I only have (non-linear) steps in frequency, not time. I am not sure as to what to provide as the sample time, but the sample time is used at multiple points in the example notebook.

Another thing I find confusing in the example is how the data is produced and passed to different function. I am unable to see how to create the right FRD object to pass to the arma function.

u = randn(10000)
I have realized, after a little back and forth, that this is data input for the simulation. This was not clear without knowledge of the lsim functionality, which I don't have. (I know this can seem like nitpicking, but I want to make aware of things that may not be clear you use as the developer, but that are apparent as a user).

Gd = c2d(Gtest, h)
y,x,t = lsim(Gd, u);

Is y the gain, x the phase and t the corresponding time? I am uncertain about things.

y is defined from a lsim simulation of the discrete version of the real system, that is fine lsim simulates in time, based on the third line of it's docstring. But then y is data from a time-simulation, right? y is then passed, in transposed form, to iddata:
d = iddata(y', u', h)
d is now also dependent on time-data, right?

This variable d is then passed to tfest
H = tfest(d, 0.05)[1]

So the new variable H is indirectly dependent on time-data, right? And H is what is passed to arma. So where is this process would one provide the frequency response data?

The code run as number 17 in https://nbviewer.jupyter.org/github/JuliaControl/ControlExamples.jl/blob/master/identification_time_vs_freq.ipynb does not run for me, and produces the following error as the final output when I pase it line-by-line into the REPL (with the packeges loaded):

julia> Gtest = tf([12], [1.1, 0.0009, 12]) # True system
TransferFunction{Continuous,ControlSystems.SisoRational{Float64}}
          12.0
-------------------------
1.1*s^2 + 0.0009*s + 12.0

Continuous-time transfer function model

julia> h = 0.01                            # Sample time
0.01

julia> u = randn(10000)
10000-element Array{Float64,1}:
  0.024656905099496182
  0.7369168037019351
 -0.9423553147187463
  0.37136329293201015
 -0.3167506517388129
 -0.0857300799373316
  0.7012325268683016
 -0.4227563751716908
  0.07708947019985367
 -0.5844135141678407
  0.4594979454511769
  0.4939273214442418
  0.9955279747485412
  ⋮
 -1.1531717187843091
  0.7344591282679885
 -0.7862105956710846
 -0.11295530888595104
 -1.1619057490283533
  0.22700580972098328
 -0.6349170652373435
 -0.4580489013680076
  0.5386217443220899
 -0.49982157631397645
 -0.7563244661321952
 -0.9583347851245462

julia> Gd = c2d(Gtest, h)
TransferFunction{Discrete{Float64},ControlSystems.SisoRational{Float64}}
  0.0005454034730425494*z + 0.0005454019855269943
---------------------------------------------------
1.0*z^2 - 1.9989010127567197*z + 0.9999918182152893

Sample Time: 0.01 (seconds)
Discrete-time transfer function model

julia> y,x,t = lsim(Gd, u);
ERROR: MethodError: no method matching lsim(::TransferFunction{Discrete{Float64},ControlSystems.SisoRational{Float64}}, ::Array{Float64,1})
Closest candidates are:
  lsim(::TransferFunction, ::Any, ::Any, ::Any...; kwargs...) at C:\Users\densb\.julia\packages\ControlSystems\fDg6j\src\timeresp.jl:178
  lsim(::Any, ::Any, ::Any, ::Any) at deprecated.jl:70
  lsim(::Any, ::Any, ::Any, ::Any, ::Any) at deprecated.jl:70
  ...
Stacktrace:
 [1] top-level scope at REPL[21]:1

So I am unable to make use of the implementation, as I don't understand some important things in using it. If things become clear to me, I would love to help in making the examples more clear.

And finally, I want to say THANK YOU. It is crazy to me how fast you simply implemented what I asked for, updated examples etc. I can't say anything about you in real life, but as a package manager on GitHub, you are a saint <3

Answer 7 · 2020-10-24T09:43:53.000Z

Everything related to the time domain in the notebook is simply there to generate some fake data, and compare the frequency domain estimation to the time domain estimation. If you already have frequency domain data, what you have is H in the notebook. This is an object of type FRD which contains a vector of frequencies and the complex frequency response, similar to in matlab. If you have magnitude and phase, you need to convert those to complex numbers
z = mag * cis(phase).

Thanks for your kind words 😃

Answer 8 · 2020-10-25T06:02:47.000Z

I realized that arma is a terrible name for this function, it does not estimate an ARMA model at all. I changed the name to tfest like I believe it´s called in matlab. tfest thus caliculates a nonparametric model if fed an IdData object, and a parametric model if fed an FRD object and an initial guess.

Answer 9 · 2020-10-25T10:22:58.000Z

and a parametric model if fed an FRD object and an initial guess.

I am trying to create a vector of frequencies, and a vector of the response as a complex number, and pass it to the contructor FRD to create an FRD object to pass to tfest along with an initial guess. But I am getting the following error

ERROR: MethodError: no method matching output(::FRD{Array{Float64,1},Array{Complex{Int64},1}})
Closest candidates are:
  output(!Matched::ControlSystemIdentification.AbstractIdData) at C:\Users\densb\.julia\packages\ControlSystemIdentification\MhFCT\src\types.jl:58
  output(!Matched::AbstractArray) at C:\Users\densb\.julia\packages\ControlSystemIdentification\MhFCT\src\types.jl:61

after running the last line in this MWE:

## Minimal working example:
using ControlSystemIdentification
frequencies = [1.; 2; 3]
responses = [1+1im; 2+2im; 3+3im]
data_as_FRD = FRD(frequencies, responses)

guess = tf([1, 2], [1, 2, 3])

tfest(data_as_FRD, guess)

I can not see why tfest does not return an estimated transfer function when is given an FRD object and an initial guess. Do you know where this goes wrong?

Answer 10 · 2020-10-25T10:44:14.000Z

Have you checked out the master branch? I haven't tagged a new release yet

Answer 11 · 2020-10-25T10:50:01.000Z

Have you checked out the master branch? I haven't tagged a new release yet

No, I have simply updated the package and tried it. So I should try to copy some code from the master branch manually? I'm going looking now ^_^

Sorry for requiring this much hand-holding, and thank you a lot for it. Once I have it working, I intend to create a notebook example for when you have the frequency response data.

And on that note - is it not better to document in Pluto-notebooks than in Jupter, as Pluto has "everything can be seen in the notebook" as a foundational principle? This means that all code needed to get the results is in the notebook, guaranteed. This is an error-finder for the documenter and a guarantee for the user. I was just thinking that all notebook documentation should be Pluto due to this extra security... thoughts?

Answer 12 · 2020-10-25T11:00:30.000Z

Just add the master branch using the package manager.

Pluto notebooks do not store and show the output, so they are not useful for documentation unless you check in an html file or pdf as well.

They also do not guarantee anything, if I call local code on my machine no one else can reproduce it.

Answer 13 · 2020-10-25T13:58:42.000Z

Just add the master branch using the package manager.

Done. I did not include "master" anywhere in the URL, but it seemed to work fine. I simply added "https://github.com/baggepinnen/ControlSystemIdentification.jl". I believe that should be the master branch, as it is shown by default when one navigates to the URL

I am now getting an inexact error from using tfest in the last line of my MWE: ERROR: InexactError: Int64(NaN). I also get this error when I am trying to enter real data, so I don't think that the problem is that the data is unrealistic. I have entered my real data for magnitude as dB, log10 and simply the actual gain as I was uncertain about which format the function takes, but those three attempts AND the MWE produce the same inexact error when attempting to call tfest (which is, as you correctly said, what matlab calls it. Nice).

I just want to add that the problem is non-critical, and I don't want you to spend your entire Sunday helping me out unless you really want to. You can check it out when you want, and I am very grateful for it whenever that is.

Pluto notebooks do not store and show the output, so they are not useful for documentation unless you check in an html file or pdf as well.

I was thinking that one would save the notebook as a static HTML only, yes. I don't know if it is common to run the Jupyter examples, but I have only ever looked at them, so I and HTML with be same-same for me.

They also do not guarantee anything, if I call local code on my machine no one else can reproduce it.

True. But if one would refrain from that, as one should in examples, it should be good...

EDIT: In case it is relevant, here is the version info:

[[ControlSystemIdentification]]
deps = ["BandedMatrices", "ComponentArrays", "ControlSystems", "DSP", "FFTW", "FillArrays", "LinearAlgebra", "LowLevelParticleFilters", "MatrixEquations", "MonteCarloMeasurements", "Optim", "Parameters", "Random", "RecipesBase", "Roots", "Statistics", "StatsBase", "TotalLeastSquares"]
git-tree-sha1 = "5ebaa2816ce7dc2e22536067d7e4b749a3508d1e"
repo-rev = "master"
repo-url = "https://github.com/baggepinnen/ControlSystemIdentification.jl"
uuid = "3abffc1c-5106-53b7-b354-a47bfc086282"
version = "1.2.1"

Answer 14 · 2020-12-12T15:18:10.000Z

I've used the new tfest for a few weeks now and it appears to be working quite well. If you are struggling with the interface, feel free to open a discussion