armax.jl 3.22 KB
Newer Older
1
2
using Convex
using SCS
3
"""`ar(y, na; λ = 0)`"""
4
function ar(y::Vector{Float64}, na; λ = 0, normtype=2, verbose=false)
5
    y_train, A = getARregressor(y,na)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
6
    n_points = size(A,1)
7
8
9
    if normtype == 2
        if λ == 0
            w = A\y_train
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
10
            method = :LS
11
12
        else
            w = (A'A + λ*eye(size(A,2)))\A'y_train
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
13
            method = :LS_reg
14
15
16
17
18
19
        end
    elseif normtype == 1
        w = Variable(size(A,2),1)
        problem = minimize(sum(abs(A*w-y_train )) + λ*norm(w))
        solve!(problem, SCSSolver(verbose=Int(verbose)))
        w = w.value[:]
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
20
        method = :L1
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
21
22
23
    end
    prediction = A*w
    error = y_train - prediction
24
    model = AR(w,na)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
25
    result = FitResult(y_train,prediction,na, method,λ)
26
    return model, result
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
27
end
28
ar(iddata::IdData, na; λ = 0) = ar(iddata.y, na; λ = 0)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
29

30
"""arx(y, u, na, nb; λ = 0)"""
31
function arx(y::Vector{Float64}, u::VecOrMat{Float64}, na, nb; λ = 0, normtype=2, verbose=false)
32
    y_train, A = getARXregressor(y,u, na, nb)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
33
    n_points = size(A,1)
34
35
36
    if normtype == 2
        if λ == 0
            w = A\y_train
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
37
            method = :LS
38
39
        else
            w = (A'A + λ*eye(size(A,2)))\A'y_train
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
40
            method = :LS_reg
41
42
43
44
45
46
        end
    elseif normtype == 1
        w = Variable(size(A,2),1)
        problem = minimize(sum(abs(A*w-y_train )) + λ*norm(w))
        solve!(problem, SCSSolver(verbose=Int(verbose)))
        w = w.value[:]
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
47
        method = :L1
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
48
49
50
    end
    prediction = A*w
    error = y_train - prediction
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
51
52

    si = na+1
53
    b = Polynom[w[si:si+nb[1]-1]]
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
54
55
    si += nb[1]
    for i = 2:length(nb)
56
        push!(b,w[si:si+nb[i]-1])
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
57
58
        si += nb[i]
    end
59
    model = ARX(w[1:na],b,na,nb)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
60
    result = FitResult(y_train,prediction,na, method, λ)
61
    return model, result
62
end
63
arx(iddata::IdData, na, nb; λ = 0) = arx(iddata.y, iddata.u, na, nb; λ = 0)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
64
65


66
function getARregressor(y::Vector{Float64},na)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
67
68
69
70
71
72
    A = toeplitz(y[na+1:end],y[na+1:-1:1])
    y = copy(A[:,1])
    A = A[:,2:end]
    n_points = size(A,1)
    return y,A
end
73
getARregressor(iddata::IdData, na) = getARXregressor(iddata.y, na)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
74

75
function getARXregressor(y::Vector{Float64},u::VecOrMat{Float64}, na, nb)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
76
    assert(length(nb) == size(u,2))
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
77
78
79
80
81
82
    m = max(na+1,maximum(nb))
    n = length(y) - m+1
    offs = m-na-1
    A = toeplitz(y[offs+na+1:n+na+offs],y[offs+na+1:-1:1])
    y = copy(A[:,1])
    A = A[:,2:end]
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
83
84

    for i = 1:length(nb)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
85
        offs = m-nb[i]
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
86
87
88
        A = [A toeplitz(u[nb[i]+offs:n+nb[i]+offs-1,i],u[nb[i]+offs:-1:1+offs,i])]
    end

89

Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
90
91
    return y,A
end
92
getARXregressor(iddata::IdData, na, nb) = getARXregressor(iddata.y,iddata.u, na, nb)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
93
94
95


"""Plots the RMSE and AIC for model orders up to `n`. Useful for model selection"""
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
96
function find_na(y::Vector,n::Int)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
97
98
99
100
    error = zeros(n,2)
    for i = 1:n
        w,e = ar(y,i)
        error[i,1] = rms(e)
101
        error[i,2] = aic(e,i)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
102
103
104
        print(i,", ")
    end
    println("Done")
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
105
    scatter(error, show=true)
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
106
end
107

Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
108
109
110
111
112
113
"""
    plotmodel(y,m::AR)

Plots a signal `y` and the output of the model `m`
"""
function plotmodel(y,m::AR)
114
115
116
117
    na = length(m.a)
    y,A = getARregressor(y,na)
    yh = A*m.a
    error = y-yh
Fredrik Bagge Carlson's avatar
Fredrik Bagge Carlson committed
118
119
120
121
    plot(y,c=:black)
    plot(yh,c=:b)
    plot(error,c=:r, title="Fitresult, AR, n_a: $na, RMSE = $(rms(error)) Fit = $(fit(y,yh))")

122
end