From e765c299ccc2a17a876911f8b8f1807309edc2d2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ram=C3=B3n=20Calvo=20Gonz=C3=A1lez?= <ramon@ramoncalvo.com>
Date: Sat, 10 Dec 2022 12:09:43 +0100
Subject: [PATCH] Exercises 7 & 8

---
 exercise-7/exercise.m | 116 ++++++++++++++++++++++++++++++++++++++++++
 exercise-8/exercise.m |  66 ++++++++++++++++++++++++
 2 files changed, 182 insertions(+)
 create mode 100644 exercise-7/exercise.m
 create mode 100644 exercise-8/exercise.m

diff --git a/exercise-7/exercise.m b/exercise-7/exercise.m
new file mode 100644
index 0000000..46f7c13
--- /dev/null
+++ b/exercise-7/exercise.m
@@ -0,0 +1,116 @@
+clear, close, clc
+
+%% Exercise 7.1
+a = [0.25, -0.2, 0.1, 0.05];
+b = [0.6, 0.3, -0.05];
+N = 100;
+sigma_u = 1;
+sigma_d = 0.1;
+
+u_id = randn(N, 1) * sqrt(sigma_u);
+[y, y_noiseless] = generate_dataset(a, b, u_id, sigma_d, N);
+
+
+%% Exercise 7.2
+
+theta_id = cell(8, 1);
+for i=1:8
+    theta_id{i} = compute_arx(y, u_id, i);
+end
+
+%% Exercise 7.3
+
+N_val = 100;
+sigma_u = 1;
+sigma_d_val = 0.1;
+
+u_val = randn(N_val, 1) * sqrt(sigma_u);
+[y_val, y_val_noiseless] = generate_dataset(a, b, u_val, sigma_d_val, N_val);
+
+%% Exercise 7.4
+
+plot(y_val_noiseless, "--black", 'LineWidth', 2);
+title('Predicted validation output')
+xlabel('Timestep')
+ylabel('y')
+hold on
+plot(y_val, "blacko")
+y_val_pred = cell(8, 1);
+y_id_pred = cell(8, 1);
+for i=1:8
+    a = theta_id{i}(1:i);
+    b = theta_id{i}(i+1:end);
+    
+    y_val_pred{i} = inference(a, b, u_val);
+    y_id_pred{i} = inference(a, b, u_id);
+
+    plot(y_val_pred{i}, 'color', rand(1,3))
+end
+
+%% Exercise 7.5
+
+err_id = zeros(8, 1);
+err_val = zeros(8, 1);
+for i=1:8
+    err_id(i) = mean((y_id_pred{i} - y_noiseless) .^ 2);
+    err_val(i) = mean((y_val_pred{i} - y_val_noiseless) .^ 2);
+end
+figure
+plot(err_id, '--redo')
+hold on
+plot(err_val, '--greeno')
+
+%% Exercise 7.6 and 7.7 I'm just too lazy
+
+function [y, y_noiseless]=generate_dataset(a, b, u, sigma_d, N)
+    y = zeros(N, 1);
+    y_noiseless = zeros(N, 1);
+
+    d_id = randn(N, 1) * sqrt(sigma_d);
+
+    for i=1:N
+        for j=1:length(a)
+            if i - j > 0
+                y(i) = y(i) + a(j) * y(i-j);
+            end
+        end
+        for j=1:length(b)
+            if i - j > 0
+                y(i) = y(i) + b(j) * u(i-j);
+            end
+        end
+        y_noiseless(i) = y(i);
+        y(i) = y(i) + d_id(i);
+    end
+end
+
+function theta=compute_arx(y, u, order)
+    % Build the phi matrix
+    if order == 1
+        phi = [y(1:end-1) u(1:end-1)];
+    else
+        phi_y = toeplitz(y(1:end-1), zeros(order, 1));
+        phi_u = toeplitz(u(1:end-1), zeros(order, 1));
+        phi = [phi_y, phi_u];
+    end
+    
+    % Identify the system
+    theta = phi \ y(2:end);
+end
+
+function y_pred=inference(a, b, u)
+    N = length(u);
+    y_pred = zeros(N, 1);
+    for i=1:N
+        for j=1:length(a)
+            if i - j > 0
+                y_pred(i) = y_pred(i) + a(j) * y_pred(i - j);
+            end
+        end
+        for j=1:length(b)
+            if i - j > 0
+                y_pred(i) = y_pred(i) + b(j) * u(i - j);
+            end
+        end
+    end
+end
\ No newline at end of file
diff --git a/exercise-8/exercise.m b/exercise-8/exercise.m
new file mode 100644
index 0000000..e41837e
--- /dev/null
+++ b/exercise-8/exercise.m
@@ -0,0 +1,66 @@
+clear, close, clc
+
+N = 10000;
+
+a = [1, -1.5, 0.7];
+b = [0.0, 1.0, 0.5];
+c = [1.0, -1.0, 0.2];
+
+B_A = tf(b, a, 1.0);
+C_A = tf(c, a, 1.0);
+
+e = randn(N, 1);
+e_u = randn(N, 1);
+
+u = generate_input_sequence(e, e_u);
+
+y = lsim(B_A, u) + lsim(C_A, u);
+
+%% Exercise 8.1
+
+C_inv = tf([1, 0, 0], c, 1.0);
+
+u_F = lsim(C_inv, u);
+y_F = lsim(C_inv, y);
+
+phi_u_F = toeplitz(u_F(1:end-1), zeros(2, 1));
+phi_y_F = toeplitz(-y_F(1:end-1), zeros(2, 1));
+phi = [phi_y_F, phi_u_F];
+
+theta = phi \ y_F(2:end)
+% G_arx = arx([y_F u_F], [2 2 1])
+
+%% Exercise 8.2
+
+N_val = 1000;
+y_val = y(end-N_val:end);
+u_val = u(end-N_val:end);
+
+G_pred = tf([theta(3), theta(4)], [1.0, theta(1), theta(2)], 1.0);
+y_pred = lsim(G_pred, u_val);
+
+plot(y_val, '--black');
+title('Validation plot')
+xlabel('Time')
+ylabel('Output')
+hold on
+plot(y_pred, 'r')
+
+%% Exercise 8.3 is too ez and I'm mega lazy
+
+%% Helper functions
+function u=generate_input_sequence(e, e_u)
+    N = size(e, 1);
+    u = zeros(N, 1);
+
+    for i=1:N
+        if i > 1
+            u(i) = u(i) + 0.1 * u(i - 1);
+        end
+        if i > 2
+            u(i) = u(i) + 0.12 * u(i - 2);
+        end
+
+        u(i) = u(i) + e(i) + e_u(i);
+    end
+end