baggepinnen / ControlSystemIdentification.jl

System Identification toolbox, compatible with ControlSystems.jl
https://baggepinnen.github.io/ControlSystemIdentification.jl/dev
MIT License
132 stars 13 forks source link

arx() estimator issues #75

Closed peroveh closed 3 years ago

peroveh commented 3 years ago

Apparent issues with several of the arx() estimators used in your identification_arx.ipynb example with slightly modified AR and ARMA transfer functions with the following solvers: (I assume they all should handle the MA part also given nb is specified? I did not find an armax version)

Gls      = arx(d,na,nb)                      # Regular least-squares estimation
Gtls     = arx(d,na,nb, estimator=tls)       # Total least-squares estimation
Gwtls    = arx(d,na,nb, estimator=wtls_estimator(y,na,nb)) # Weighted Total least-squares estimation
Gplr, Gn = plr(d,na,nb,nc, initial_order=20) # Pseudo-linear regression
Gn4sid = tf(n4sid(d, :auto; verbose=false).sys) #

1) Tested pure AR with two imaginary poles on

     0.8
-------------
1.0z^2 + 0.25

A small error is added to simulated data

y  = sim(G,u)

eps = 0.000001
yn = y + eps * randn(size(y));

Testing various specifications of system order (denominator and numerator). I am a bit uncertain if n is the number of coefficients here or the order?

With na=2,nb=2,nc=1 (?) they all gave reasonably similar and ok results..

Gwtls = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
     -1.830761022395132e-8z + 0.8000000364406703
-----------------------------------------------------
1.0z^2 + 1.1336180899546022e-7z + 0.24999999923772648

Gn4sid = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
1.0482460250016439e-7z^2 - 2.4016742360023076e-8z + 0.8000000313485508
----------------------------------------------------------------------
          1.0z^2 + 7.743629093504012e-8z + 0.250000013686116

With na=2, nb=1 (should suffice for the numerator) the arx versions start failing bad, while n4sid still does well

Gls = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
              -0.008425068831949092
-------------------------------------------------
1.0z^2 - 0.0149157451820707z + 0.2485242562797274

Gtls = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
               -0.14710343635360346
--------------------------------------------------
1.0z^2 - 0.10034853554694595z + 1.0435201347058614

Gplr = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
               -0.00843907141354728
---------------------------------------------------
1.0z^2 - 0.01502858249704049z + 0.25577236172028245

Gwtls = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
               -0.06823555574414619
--------------------------------------------------
1.0z^2 - 0.09093792215356927z + 0.9887132202347376

Gn4sid = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
-7.035087559097725e-8z^2 + 1.8586494326913326e-8z + 0.7999999720915703
----------------------------------------------------------------------
         1.0z^2 - 4.118840077338426e-8z + 0.2500000350249592

With na=nb=4 (too high?) , the arx and plr also goes bad (large errs in both numerator and denom) Again n4sid good.

Gls = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
     4.498197651604916e-8z^3 + 0.7999999908664831z^2 + 0.003570318828339773z - 0.11531636330533376     
-------------------------------------------------------------------------------------------------------
1.0z^4 + 0.0044629928399430585z^3 + 0.10585461073330074z^2 + 0.001115779335659886z - 0.0360363583918387

Gtls = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
    2.7838920503841905e-8z^3 + 0.8000000588632808z^2 + 0.891535197737989z + 0.6889073529493018
--------------------------------------------------------------------------------------------------     
1.0z^4 + 1.1144190526902322z^3 + 1.111134362289782z^2 + 0.27860493890406324z + 0.21528368611087234     

Gplr = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
      4.420389876275559e-8z^3 + 0.8000000013003751z^2 - 0.003497016602469035z - 0.1134473493326818
---------------------------------------------------------------------------------------------------------
1.0z^4 - 0.0043711788328240336z^3 + 0.10819087189424594z^2 - 0.0010927728758850003z - 0.03545228883099305

Gwtls = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
  1.5679542424341932e-8z^3 + 0.8000000093945829z^2 + 0.038695551344337255z + 0.45693930399129623
--------------------------------------------------------------------------------------------------
1.0z^4 + 0.04836949844579501z^3 + 0.821174204714606z^2 + 0.01209246301841781z + 0.1427935848734602

Gn4sid = TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Float64}}
1.8043888767223154e-8z^2 + 6.290335821529143e-8z + 0.8000000000517333
---------------------------------------------------------------------
        1.0z^2 + 7.375779662033288e-8z + 0.25000005919139123

As can be seen the n4sid has typical coefficient errors in the e-8 range, while the others are up to e-1 Per Ove Husøy

baggepinnen commented 3 years ago

@mapi1 do you think this is related to https://github.com/baggepinnen/ControlSystemIdentification.jl/commit/890a3ee9e022055a136c1128146846eb402b6e14 or https://github.com/baggepinnen/ControlSystemIdentification.jl/commit/2d212e3106670708092bd4bf50c09cdc6314db8d ?

Here's some code to roughly reproduce what's seen above since no code was posted

using ControlSystemIdentification, ControlSystems, Plots
using Random, LinearAlgebra
default(size=(500,280))
N  = 500      # Number of time steps
t  = 1:N
Δt = 1        # Sample time
u  = randn(N) # A random control input
G  = tf(0.8, [1,0,0.25], 1) # An interesting system
sim(sys,u) = lsim(sys, u, t)[1][:]
y  = sim(G,u)
yn = y + randn(size(y));

# Validation data
uv  = randn(N)
yv  = sim(G,uv)
ynv = yv + randn(size(yv));

na,nb = 2,1   # Number of polynomial coefficients to estimate
d = iddata(yn,u,Δt)
Gls = arx(d,na,nb)

Chaning the inputdelay parameter does recover a better transfer function (but does not change the numerator order as would be expected)

julia> Gls = arx(d,2,1)
TransferFunction{Discrete{Int64}, SisoRational{Float64}}
               -0.053195176045275815
----------------------------------------------------
1.0z^2 + 0.020944097150896568z + 0.13061559216712304

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

julia> Gls = arx(d,2,1, inputdelay=2)
TransferFunction{Discrete{Int64}, SisoRational{Float64}}
                0.8711130327440996
---------------------------------------------------
1.0z^2 - 0.008279348253820712z + 0.1078308656207146

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

@peroveh note that arx is expected to perform poorly since what you are simulating is an output-error process but you're estimating an input-error model. They should not perform this poorly though, and I believe it's related to an internal index shift that has gone wrong.

peroveh commented 3 years ago

I was surprised that the error was low if na=nb=2 in the solver, but worse when na=2 and nb=1 (which should be sufficient)

Btw: is n the order of the largest z term, or is it the number of coefficients...? I.e. for numerator 1.1z^2+0.8z+0.1 should this require nb=2 or 3?

Den 14. august 2021 kl. 16.42.25 +02.00 skrev Fredrik Bagge Carlson @.***>:

@mapi1 https://github.com/mapi1 do you think this is related to 890a3ee

or 2d212e3 ?

Here's some code to roughly reproduce what's seen above since no code was posted

using ControlSystemIdentification, ControlSystems, Plots

using Random, LinearAlgebra

default(size=(500,280))

N = 500 # Number of time steps

t = 1:N

Δt = 1 # Sample time

u = randn(N) # A random control input

G = tf(0.8, [1,0,0.25], 1) # An interesting system

sim(sys,u) = lsim(sys, u, t)[1][:]

y = sim(G,u)

yn = y + randn(size(y));

Validation data

uv = randn(N)

yv = sim(G,uv)

ynv = yv + randn(size(yv));

na,nb = 2,1 # Number of polynomial coefficients to estimate

d = iddata(yn,u,Δt)

Gls = arx(d,na,nb) Chaning the inputdelay parameter does recover a better transfer function julia> Gls = arx(d,2,1)

TransferFunction{Discrete{Int64}, SisoRational{Float64}}

-0.053195176045275815


1.0z^2 + 0.020944097150896568z + 0.13061559216712304

Sample Time: 1 (seconds)

Discrete-time transfer function model

julia> Gls = arx(d,2,1, inputdelay=2)

TransferFunction{Discrete{Int64}, SisoRational{Float64}}

0.8711130327440996


1.0z^2 - 0.008279348253820712z + 0.1078308656207146

Sample Time: 1 (seconds)

Discrete-time transfer function model @peroveh https://github.com/peroveh note that arx is expected to perform poorly since what you are simulating is an output-error process but you're estimating an input-error model. They should not perform this poorly though, and I believe it's related to an internal index shift that has gone wrong.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/baggepinnen/ControlSystemIdentification.jl/issues/75#issuecomment-898902872, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD7ALMGX2X6SVTBEEVMGWCLT4Z6FDANCNFSM5CAIMU5A. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email.

mapi1 commented 3 years ago

@mapi1 do you think this is related to 890a3ee or 2d212e3 ?

I don't think it is particularly related, but it probably made things more complicated.

I also think this behavior is correct (though probably not well enough documented). To clarify this a little bit: nb is the number of coefficients, while the order is given by nb + inputdelay - 1, which also differs from the definition used by eg. Matlab potentially confusing users coming from there. So for the example above:

G = tf(0.8, [1,0,0.25], 1)

TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Float64}}
0.8
-------------
1.0z^2 + 0.25

nb = 1 is enough with respect to the number of coefficients but as the order of B is 2, the estimation still fails. Setting inputdelay = 2 fixes this while maintaining the number of coefficients.

Maybe it would help to print something like the following in order to make clear what is actually estimated?

julia> Gls = arx(d,2,1)
┌ Info: Estimating TransferFunction of form: 
│        b₁z^1 
│  ------------------- 
└  1.0z^2 + a₁z^1 + a₂
...

julia> Gls = arx(d,2,1,inputdelay = 2)
┌ Info: Estimating TransferFunction of form: 
│          b₁
│  ------------------- 
└  1.0z^2 + a₁z^1 + a₂
...
baggepinnen commented 3 years ago

There's probably still an error though, as the z^1 in your printout is not reflected in the actual transfer functions returned by the estimation

julia> Gls = arx(d,2,1)
               -0.053195176045275815
----------------------------------------------------
1.0z^2 + 0.020944097150896568z + 0.13061559216712304

julia> Gls = arx(d,2,1, inputdelay=2)
                0.8711130327440996
---------------------------------------------------
1.0z^2 - 0.008279348253820712z + 0.1078308656207146

no matter how the inputdelay is defined, there should be a difference in the order of the z in the numerator here?

mapi1 commented 3 years ago

You are right, there is an error. The case na > nb, (or better order of A > order of B) is missing in params2poly2 and from the tests. I will open a PR fixing this and trying to be more clear in the documentation.

Apart from that, I found something that is strange, this little snippet passes as true when I add it to the tests, though it shouldn't right? What am I missing?

G1  = tf([0.8], [1,0,0.25], 
G2  = tf([0.8, 0], [1,0,0.25], 1) 
@test isapprox(G1, G2)
baggepinnen commented 3 years ago

I will open a PR fixing this and trying to be more clear in the documentation.

Very much appreciated!

though it shouldn't right?

You're right, that appears to be a bug in ControlSystems.jl

baggepinnen commented 3 years ago

Thanks to @mapi1, this issue is now resolved on the master branch of this package. To get your desired behavior, you must supply the correct value to the inputdelay keyword argument, e.g.,

julia> Gls = arx(d,na,nb, inputdelay=2)
TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Float64}}
                 0.8368547377892526
----------------------------------------------------
1.0z^2 - 0.002375293735401682z + 0.07787696411618848
peroveh commented 3 years ago

Hi Fredrik I am a bit uncertain about these inputs.. Given that na is the (maximum?) order of the denominator, nb is of the numerator,,, what is the extra inputdelay? Will this be the max of na and nb? Or is the model a concatenation of a arma and a fixed delay of two timesteps? Are there any docomentation to look at for usage and algorithmic description?

Per Oe On August 20, 2021 at 2:15:17 pm +02:00, Fredrik Bagge Carlson @.***> wrote:

Thanks to @mapi1 https://github.com/mapi1, this issue is now resolved on the master branch of this package. To get your desired behavior, you must supply the correct value to the inputdelay keyword argument, e.g.,

julia> Gls = arx(d,na,nb, inputdelay=2) TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Float64}} 0.8368547377892526

1.0z^2 - 0.002375293735401682z + 0.07787696411618848

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/baggepinnen/ControlSystemIdentification.jl/issues/75#issuecomment-902649637, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD7ALMEOPQP4V7QLDTXKNE3T5ZBNLANCNFSM5CAIMU5A. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email.

peroveh commented 3 years ago

Thanks to @mapi1, this issue is now resolved on the master branch of this package. To get your desired behavior, you must supply the correct value to the inputdelay keyword argument, e.g., How do I ensure that i get the corrected vesion in the master branch? Can I specify a specific version of ControlSystems.jl or do I need to download the latest Julia release?

mapi1 commented 3 years ago

To install the latest version you can just use ] + add ControlSystemIdentification#master from the REPL.

In the latest version, the documentation help?> arx is then also updated which should answer some of your questions.

baggepinnen commented 3 years ago

Using the master branch will also require using the master branch of ControlSystems, as there has been some yet-to-be-released breaking changes.