From 39f96b244b1bcca46b963ee6969935e2839fc5af Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <fredrikb@control.lth.se>
Date: Sun, 6 Sep 2015 19:23:43 +0200
Subject: [PATCH] Implemented minimization of 1-norm for ar arx
---
src/PCA.jl | 9 ++-------
src/armax.jl | 36 ++++++++++++++++++++++++++----------
2 files changed, 28 insertions(+), 17 deletions(-)
diff --git a/src/PCA.jl b/src/PCA.jl
index dc290a9..c852e7d 100644
--- a/src/PCA.jl
+++ b/src/PCA.jl
@@ -1,13 +1,8 @@
-""" Performs PCA PCA(W,doplot=false)"""
-function PCA(W,doplot=false)
+""" Performs PCA PCA(W)"""
+function PCA(W)
W0 = mean(W,1);
W = W-repmat(W0,size(W,1),1);
(score,latent,C) = svd(W)
score = score*diagm(latent)
- if doplot
- newplot(cumsum(latent)./sum(latent),"o")
- title("Fraction of variance explained")
- ylim([0,1])
- end
C,score,latent,W0
end
diff --git a/src/armax.jl b/src/armax.jl
index 7dc109f..9773a26 100644
--- a/src/armax.jl
+++ b/src/armax.jl
@@ -1,11 +1,20 @@
+using Convex
+using SCS
"""`ar(y, na; λ = 0)`"""
-function ar(y::Vector{Float64}, na; λ = 0)
+function ar(y::Vector{Float64}, na; λ = 0, normtype=2, verbose=false)
y_train, A = getARregressor(y,na)
n_points = size(A,1)
- if λ == 0
- w = A\y_train
- else
- w = (A'A + λ*eye(size(A,2)))\A'y_train
+ if normtype == 2
+ if λ == 0
+ w = A\y_train
+ else
+ w = (A'A + λ*eye(size(A,2)))\A'y_train
+ 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[:]
end
prediction = A*w
error = y_train - prediction
@@ -16,13 +25,20 @@ end
ar(iddata::IdData, na; λ = 0) = ar(iddata.y, na; λ = 0)
"""arx(y, u, na, nb; λ = 0)"""
-function arx(y::Vector{Float64}, u::VecOrMat{Float64}, na, nb; λ = 0)
+function arx(y::Vector{Float64}, u::VecOrMat{Float64}, na, nb; λ = 0, normtype=2, verbose=false)
y_train, A = getARXregressor(y,u, na, nb)
n_points = size(A,1)
- if λ == 0
- w = A\y_train
- else
- w = (A'A + λ*eye(size(A,2)))\A'y_train
+ if normtype == 2
+ if λ == 0
+ w = A\y_train
+ else
+ w = (A'A + λ*eye(size(A,2)))\A'y_train
+ 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[:]
end
prediction = A*w
error = y_train - prediction
--
GitLab