diff --git a/code/week02/Ex1_yhofmann.m b/code/week02/Ex1_yhofmann.m new file mode 100644 index 0000000..e2f34bf --- /dev/null +++ b/code/week02/Ex1_yhofmann.m @@ -0,0 +1,46 @@ +v = randn(100,1)*sqrt(0.3) + 0.5; +f = @(theta) -sum((y - (theta(1)*x + theta(2)*(2*x.^2 -1) + v)).^2) + +theta_initial = [0,0] +options = optimset('Display', 'iter', 'TolFun', 1e-6); +estimated_theta = fminunc(f, theta_initial, options); + +fprintf('Estimated theta1: %.4f\n', estimated_theta(1)); +fprintf('Estimated theta2: %.4f\n', estimated_theta(2)); +%% + +% Number 1 + +for i = 1:size(x,1) +w_k = [x(i), 2*x(i)^(2)-1]; +w(i,1) = w_k(1); +w(i,2) = w_k(2); +end + + +%better +w = [x 2*x.^2-1]; + + + + +theta = inv(transpose(w)*w)*transpose(w)*(y-0.5) + + +% Number 2 +var_v = 0.3; +mu_v = 0.5; +var_theta =0.02; +mu_theta = [1.3; 0.9]; +theta_MAP = inv((1/var_v)*transpose(w)*w + (1/var_theta)*eye(2)) * ((1/var_v)*transpose(w)*(y-mu_v)+mu_theta/var_theta) + +%Number 3 +error_ML = 0; +error_MAP =0; +for i = 1:size(x_v,1) +error_ML = error_ML + (theta(1)*x_v(i) + theta(2)*(2*(x_v(i))^2-1) + 0.5 - y_v(i))^2; +error_MAP = error_MAP + (theta_MAP(1)*x_v(i) + theta_MAP(2)*(2*(x_v(i))^2-1) + 0.5 - y_v(i))^2; + +end +error_ML =error_ML/size(x_v,1) +error_MAP = error_MAP/size(x_v,1) \ No newline at end of file diff --git a/code/week02/SystemID_HS2023_Solution1.pdf b/code/week02/SystemID_HS2023_Solution1.pdf new file mode 100644 index 0000000..3408df9 Binary files /dev/null and b/code/week02/SystemID_HS2023_Solution1.pdf differ diff --git a/code/week03/Ex2_yhofmann.m b/code/week03/Ex2_yhofmann.m new file mode 100644 index 0000000..291aca3 --- /dev/null +++ b/code/week03/Ex2_yhofmann.m @@ -0,0 +1,110 @@ + +clear all +close all +clc + + + +%% Ex 3 +close all + +data1 = importdata("experiment1.dat").data; +t1 = data1(:,1), +y1 = data1(:,2) +std1= data1(:,3) +reg = [ones(length(y1),1), t1, t1.^2 ]; +theta = reg\y1 + + +R = diag(std1.^2); +Z_t_direct = theta /y1; % --> Not good don't do this as probably some numerical stability issues +% or with formula +Z_t_formula = inv(reg'*inv(R)*reg)*reg'*inv(R); + + + +% dont do cov_theta1 = Z_t_direct*R*transpose(Z_t_direct) + +% Both work well +cov_theta2 = Z_t_formula*R*Z_t_formula' +%or direct formula +cov_theta3 =inv(reg'*inv(R)*reg); + +cov1 = sqrt(cov_theta3(1,1)) +cov2 = sqrt(cov_theta3(2,2)) +cov3 = 2*sqrt(cov_theta3(3,3)) + +figure; +errorbar(t1,y1,std1, 'LineStyle','none','Marker','.') +hold on +t_cont = t1(1):0.01:t1(end); + +y_pred= reg*theta; +y0 = theta(1,1) +v0 = theta(2,1) +a = theta(3,1) +plot(t_cont, polyval([a v0 y0],t_cont)) + +g_exp1 = 2*a + +hold off + +%% +%%Ex 3 part 2 +data2 = importdata("experiment2.dat").data; +t2 = data2(:,1), +y2 = data2(:,2) +std2= data2(:,3) +reg2 = [ones(length(y2),1), t2, t2.^2 ]; +%w = cov with formula from ex1 +R2 = diag(std2.^2); +cov =inv(reg2'*inv(R2)*reg2); +cov12 = sqrt(cov(1,1)) +cov22 = sqrt(cov(2,2)) +cov32 = 2*sqrt(cov(3,3)) + +w = R2; + + +%or for weighted: +theta2 = inv(reg2'*inv(w)*reg2)*reg2'*inv(w)*y2 +% + + + +figure; +errorbar(t2,y2,std2, 'LineStyle','none','Marker','.') +hold on + + + +y02 = theta2(1,1) +v02 = theta2(2,1) +a2 = theta2(3,1) +t_cont = t2(1):0.01:t2(end); +plot(t_cont, polyval([a2 v02 y02],t_cont)) + +g_exp2 = 2*a2 + +hold off + +%% Ex 4 +w0 = 10; +theta0 = pi/8; +data3 = importdata("experiment3.dat").data; +t3 = data3(:,1); +y3 = data3(:,2); + +figure; +plot(t3,y3) + +std3= data2(:,3); +reg3 = [ones(size(t3,1),1), t3, 0.5* t3.^2, cos(theta0 + w0*t3)]; +theta3 = reg3\y3; + +disp(theta3') + +true_val = [12.42, 44.19, -6.42, 1.53]; +disp("True values") +disp(true_val) + diff --git a/code/week03/SystemID_HS2023_Solution2.pdf b/code/week03/SystemID_HS2023_Solution2.pdf new file mode 100644 index 0000000..5ce1352 Binary files /dev/null and b/code/week03/SystemID_HS2023_Solution2.pdf differ diff --git a/code/week04/Ex3_yhofmann.m b/code/week04/Ex3_yhofmann.m new file mode 100644 index 0000000..bb3185e --- /dev/null +++ b/code/week04/Ex3_yhofmann.m @@ -0,0 +1,67 @@ +close all +clear all +clc +%% Test 1 +test_Signal = randn(1,100); +res_fft = fft(test_Signal,100); +res_self= fourierTransformEx3(test_Signal,100); +if(res_self ~= res_fft) + disp('False') +else + disp('Is correct') +end + +dif = sum(abs(res_self-res_fft)) +%is correct now + +%% +% a) +e = randn(1,1024); +mean(e); +var(e); + +%% +% b) +res = (abs(fft(e)).^2)/1024; + +%[pxx,f] = periodogram(e); Do not use the built in function +%freq = linspace(0,pi,1024); Don't use linspace +freq = 2*pi/1024*(0:1024-1); %Use this +idx = find(freq>=0 & freq =0 & omega< pi); + +% Uncomment to see plots of signals +%{ +figure(1) +plot(u) +hold on +plot(e) +hold on +plot(y_experimet) +legend('u', 'e', 'Y') + + +figure(2) +plot(rese) +hold on +plot(resu) +hold on +plot(y_experimet) +hold on +legend('H(e)', 'G(u)', 'Y') +%} + +%% +%b) +U = fft(u); +Y_exp = fft(y_experimet); + +% if one of them has many which are 0 is probably due to periodic input --> +% only calculate at these + + % figure(1) + % loglog(omega(idx), abs(U(idx))); +% title('U') +% figure(2) +% loglog(omega(idx), abs(Y_exp(idx))); +% title('Y_exp') + +G_est = Y_exp ./transpose(U); %Watch out as U is complex don't do U' but use transpose(U)!!!!!!!!! +G_fresp_true = squeeze(freqresp(G,omega)); + + +figure(3) +loglog(omega(idx),abs(G_est(idx))) +hold on +loglog(omega(idx),abs(G_fresp_true(idx))) +hold on +Err = G_est - G_fresp_true; +loglog(omega(idx),abs(Err(idx))) +legend('Estimate', 'True', 'Error'); +title('Magnitude') + + + + +%% +%c) +N_c = N/4; +omega_c = (2*pi/(T*N_c)*(0:N_c-1))'; +idx_c = find(omega_c>=0 & omega_c<= pi); + +% Does the same but not really good to do it like that +Y_experimet1 = fft(y_experimet(1:end/4)); +Y_experimet2 = fft(y_experimet(end/4+1:end/2)); +Y_experimet3 = fft(y_experimet(end/2+1:3*end/4)); +Y_experimet4 = fft(y_experimet(3*end/4+1:end)); +U1 = fft(u(1:end/4)); +U2 = fft(u(end/4+1:end/2)); +U3 = fft(u(end/2+1:3*end/4)); +U4 = fft(u(3*end/4+1:end)); +G1_est = Y_experimet1./transpose(U1); +G2_est = Y_experimet2./transpose(U2); +G3_est =Y_experimet3./transpose(U3); +G4_est =Y_experimet4./transpose(U4); + +%Better like that +u_res = reshape(u,N_c,4); +y_res = reshape(y_experimet,N_c,4); + +U_res = fft(u_res,[],1); +Y_res = fft(y_res,[],1); + +G_est_res = Y_res./U_res; + + +G_est_average = zeros(length(G1_est),1); +for k =1:length(G1_est) + G_est_average(k) = mean([G1_est(k);G2_est(k);G3_est(k); G4_est(k)]); +end +G_sum = G1_est+G2_est+G3_est+G4_est; +G_est_average2 = (G_sum)/4; % both ways work now + + +G_est_average2 = mean(G_est_res,2); + + +figure(5) +loglog(omega_c(idx_c),abs(G_est_average2(idx_c))) +hold on +loglog(omega(idx),abs(G_est(idx))) +hold on +loglog(omega(idx),abs(G_fresp_true(idx))) +legend('Averaged estimate', 'Estimate', 'True System'); +title('Magnitude from averaged') diff --git a/code/week05/SystemID_HS2023_Exercise4.pdf b/code/week05/SystemID_HS2023_Exercise4.pdf new file mode 100644 index 0000000..3084518 Binary files /dev/null and b/code/week05/SystemID_HS2023_Exercise4.pdf differ diff --git a/code/week06/HS2023_SysID_Exercise_05_20914032_yhofmann.m b/code/week06/HS2023_SysID_Exercise_05_20914032_yhofmann.m new file mode 100644 index 0000000..3668ee8 --- /dev/null +++ b/code/week06/HS2023_SysID_Exercise_05_20914032_yhofmann.m @@ -0,0 +1,351 @@ +function [u, Ru, y, Ruy, Ry, etfe_raw, spec_raw, etfe_smooth, spec_smooth] = HS2023_SysID_Exercise_05_20914032() + +%% Instructions +% Student submissions should be contained in this base script. +% Change the filename and function-name `12345678' to your personal Legi +% number before execution and submission. +% +%% Dependencies +% The following TA-provided files are included in the assignment will be +% required in order to form solutions: +% +% HS2023_SysID_Exercise_05_GenerateData.p +% +%% Outputs: +% u: The input (chirp) signal applied to the system u(k) +% Ru: The autocorrelation of the input (input_u) +% y: Output y(k) of the system +% Ruy: Cross correlation between u(k) and y(k) +% Ry: Autocorrelation of y(k) +% etfe_raw: Empirical Transfer Function Estimate +% spec_raw: Spectral estimate +% etfe_smooth: Smoothed ETFE by Hann window and inverse variance weighting +% spec_smooth: Smoothed spectral estimate by Hann window and inverse variance weighting +%% Policies +% +% The only allowable external toolbox to MATLAB is the Control Systems +% Toolbox. MATLAB default commands such as 'fft' or 'fftshift' are allowed. +% +% It is impermissible to other use MATLAB toolboxes, such as the Computer Vision +% Toolbox, the Econometrics Toolbox, the Deep Learning Toolbox, etc. +% +% All student-defined functions should be written at the bottom of this +% file. The homework code-submission will be this single file. +% +% The code of academic conduct remains in effect when preparing this +% homework. +% +%% Output format specification: +% The output will be stored in the class 'out'. Refer to output_ex_5.m for +% more detail about desired fields and structures + + +%% Begin routines +% Extract Legi from Filename +name=mfilename; +LegiNumber= name(end-7:end); + + + +% In midterm for plots specify:!!!!!!!!!!!!!!!!! +% Title +% Axis +% Legend + +%% problem parameters + +%time horizon of the data +T = 1024; + +%% responses to problem areas + + + +%% Part 1: Input design + +%create functions that will generate inputs under the specified +%requirements +f0 = 0; +f1 = 0.3; +k = (0:T-1); + +%chirp input +u = cos(2*pi*(f0*k+(f1-f0)*k.^2)/(2*T))'; +figure(10) + +%% Part 2: Input and Autocorrelation + + +figure(1) +[Ru, valueLag] = intcor(u,u); +stem(valueLag,Ru) % better also use plot +title('Autocorrelation of u') +% from solutions +Ru2 = (ifft(fft(u).*conj(fft(u))))/T; % To shift it to 0 use fftshift +figure(100) +plot(Ru2) + +%% Part 3: Gather output + +%execute the black-box system with computed inputs to generate noisy +%system responses +y = HS2023_SysID_Exercise_05_GenerateData(20914032, u); + + + +%% Part 4: Output Correlations + +%cross correlations between input and output +[Ruy, ~] = intcor(u,y); +% Ryu = conj(Ruy); +% from solution +Ruy2 = ifft(fft(u).*conj(fft(y)))/T; + +Ryu = intcor(y, u); + + +figure(11) +plot(Ruy) +figure(12) +plot(Ruy2) + +%autocorrelations of input +[Ry, ~] = intcor(y,y); +Ry2 = ifft(fft(y).*conj(fft(y)))/T; + + + +%% Part 5: Transfer Function Estimates + +%compute empirical transfer function estimates +etfe_raw = fft(y)./fft(u); + + +spec_raw = fft(Ryu)./fft(Ru); % WATCH OUT is yu and not uy + +omega = 2*pi/T *(0:T/2-1); +figure(2) +loglog(omega,abs(spec_raw(1:length(omega),:))); +title('Magnitude of estimated TF from spectral function') +ylim([0.1 10]) +xlim([6*10^-3 3]) + +figure(3) +semilogx(omega,angle(spec_raw(1:length(omega),:))); +title('Phase of estimated TF from spectral function') +ylim([-4 4]) + + +figure(4) +loglog(omega,abs(etfe_raw(1:length(omega),:))); +title('Magnitude of estimated TF from etfe') +ylim([0.1 10]) +xlim([6*10^-3 3]) + +figure(5) +semilogx(omega,angle(etfe_raw(1:length(omega),:))); +title('Phase of estimated TF from etfe') +ylim([-4 4]) + + +%% Part 6: Smoothed Transfer Functions + +%smooth the ETFE using windowing and inverse variance weighting +U = fft(u); +[omega6, Wg] = WfHann(20,T); % better gamma = 100 +zidx = find(omega6 ==0); +omega6 = [omega6(zidx:T); omega6(1:zidx-1)]; +Wg = [transpose(Wg(zidx:T));transpose(Wg(1:zidx-1))]; +a = U.*conj(U); +etfe_smooth = 0 *etfe_raw; + +for wn = 1:T + Wnorm = 0; + for xi = 1:T + widx = mod(xi-wn,T)+1; + etfe_smooth(wn) = etfe_smooth(wn)+Wg(widx)*etfe_raw(xi)*a(xi); + Wnorm = Wnorm +Wg(widx)*a(xi); + end + etfe_smooth(wn) = etfe_smooth(wn)/Wnorm; +end + +idx = find(omega6>0 & omega6M) + R(i) = R(i) + u(k) * y(k - h(i)-M); + + else + R(i) = R(i) + u(k) * y(k - h(i)); + end + end + end + end + R = R/M; +end + +function [omega,WHann] = WfHann(gamma,N) +%------------------------------------------------------------------ +% +% [omega,WHann] = WfHann(gamma,N) +% +% Create a frequency domain Hann window with width parameter gamma +% and data length N. The Hann window is a raised cosine. +% +% Roy Smith, 18 October, 2017. +% +% 6 November, 2017. Fixed bug in N even indexing. +% +%------------------------------------------------------------------ + +if nargin == 0, + disp('Syntax: [omega,W] = WfHann(gamma,N)') + return +elseif nargin ~= 2, + error('incorrect number of input arguments (2 expected)') + return +end + +% basic parameter checking +if length(gamma) > 1, + error('Width parameter, gamma, must be a scalar'); +end +if round(gamma) ~= gamma, + error('Width parameter, gamma, must be an integer'); +end +if gamma < 1, + error('Width parameter, gamma, must be positive'); +end +if length(N) > 1, + error('Calculation length, N, must be a scalar'); +end +if round(N) ~= N, + error('Calculation length, N, must be an integer'); +end +if N < 1, + error('Calculation length, N, must be positive'); +end + +% The simplest approach is to define the window in the time domain and +% then transform it to the frequency domain. + +lags = [floor(-N/2+1):floor(N/2)]'; +wHann = 0*lags; +idx = find(abs(lags) <= gamma); +wHann(idx) = 0.5*(1+cos(pi*lags(idx)/gamma)); + +% +zidx = find(lags==0); % index of the zero point. + +wH_raw = fft([wHann(zidx:N);wHann(1:zidx-1)]); +WHann(zidx:N) = wH_raw(1:N-zidx+1); % shift +ve freq to end +WHann(1:zidx-1) = wH_raw(N-zidx+2:N);% shift > pi freq to beginning +WHann = real(WHann); % should actually be real +omega = 2*pi/N*lags; + +return + +end +%------------------------------------------------------------------ + +function [lags,wHann] = WtHann(gamma,N) +%------------------------------------------------------------------ +% +% [lags,wHann] = WtHann(gamma,N) +% +% Create a Hann window with width parameter gamma and data length N. +% The Hann window is a raised cosine. +% +% Roy Smith, 18 October, 2017. +% +%------------------------------------------------------------------ + +if nargin == 0, + disp('Syntax: [lags,w] = WtHann(gamma,N)') + return +elseif nargin ~= 2, + error('incorrect number of input arguments (2 expected)') + return +end + +% basic parameter checking +if length(gamma) > 1, + error('Width parameter, gamma, must be a scalar'); +end +if round(gamma) ~= gamma, + error('Width parameter, gamma, must be an integer'); +end +if gamma < 1, + error('Width parameter, gamma, must be positive'); +end +if length(N) > 1, + error('Calculation length, N, must be a scalar'); +end +if round(N) ~= N, + error('Calculation length, N, must be an integer'); +end +if N < 1, + error('Calculation length, N, must be positive'); +end + +lags = [floor(-N/2+1):floor(N/2)]'; +wHann = 0*lags; +idx = find(abs(lags) <= gamma); +wHann(idx) = 0.5*(1+cos(pi*lags(idx)/gamma)); + +return + +end +%------------------------------------------------------------------ + + diff --git a/code/week06/SystemID_HS2023_Exercise5_solution.pdf b/code/week06/SystemID_HS2023_Exercise5_solution.pdf new file mode 100644 index 0000000..682f3dc Binary files /dev/null and b/code/week06/SystemID_HS2023_Exercise5_solution.pdf differ diff --git a/code/week08/Ex6_yhofmann.m b/code/week08/Ex6_yhofmann.m new file mode 100644 index 0000000..6cddc03 --- /dev/null +++ b/code/week08/Ex6_yhofmann.m @@ -0,0 +1,187 @@ +clear all +close all +%%1) +% Define the model parameters +a = 5/8; +b = 11/10; +c = 5/6; +covariance = 0.05; +G = tf([b, c], [1, -a, 0], 1); + +% Define the number of periods and time steps +num_periods = 2; +tau_max = 80; +num_samples = 100; %Magic :) +num_sinusiods = tau_max/2; + +% Initialize the input signal u(k) +u = generateSinusoidalInput(num_samples, num_periods, num_sinusiods); + + +% No noise + y_no_noise = generateOutput(u,zeros(num_samples*num_periods,1),a,b,c,num_periods,num_samples,G); + + g0 = timeEstimation(u,tau_max,y_no_noise); + + g = zeros(tau_max+1,3); + +% Simulate the system and estimate the impulse response +for realization = 1:3 + % Generate random noise + v = sqrt(covariance) * randn(num_samples*num_periods,1); + + y = generateOutput(u,v,a,b,c,num_periods,num_samples,G); + + g(:,realization) = timeEstimation(u,tau_max,y); +end + +% Calculate Average +g_hat_final = mean(g,2); + +error1 = vecnorm(g_hat_final - g0,2,1); + + + +%% +% 2) +max = 80; +g = zeros(tau_max+1,max-1); +for cur_period = 2:max + u = generateSinusoidalInput(num_samples, cur_period, num_sinusiods); + + v = sqrt(covariance) * randn(num_samples*cur_period,1); + + y = generateOutput(u,v,a,b,c,cur_period,num_samples,G); + + g(:,cur_period-1) = timeEstimation(u,tau_max,y); +end + +error = vecnorm(g - g0,2,1); +figure(1) +loglog(2:max,error); +title('Error for part 2') + +%% +% 3) +g_rand = zeros(tau_max+1,max-1); +for cur_period = 2:max + u = 4*(rand(cur_period*num_samples,1) - 0.5); + + v = sqrt(covariance) * randn(num_samples*cur_period,1); + + y = generateOutput(u,v,a,b,c,cur_period,num_samples,G); + + g_rand(:,cur_period-1) = timeEstimation(u,tau_max,y); +end + + + + +error_rand = vecnorm(g_rand - g0,2,1); +figure(2) +plot(2:max,error_rand); +hold on +plot(2:max,error) +legend('Random', 'Sinusoidal') +title('Error for part 3') + + + + +%% +% 4) +u = generateSinusoidalInput(num_samples, num_periods, num_sinusiods); +g4 = zeros(tau_max+1,200); +for realization = 1:200 + % Generate random noise + v = sqrt(covariance) * randn(num_samples*num_periods,1); + + y = generateOutput(u,v,a,b,c,num_periods,num_samples,G); + + if realization > 1 + g4(:,realization) = (timeEstimation(u,tau_max,y) + g4(:,realization-1)); + else + g4(:,realization) = timeEstimation(u,tau_max,y); + end +end + +for i = 1:size(g4,2) + g4(:,i) = g4(:,i)/i; +end + +error4 = vecnorm(g4 - g0,2,1); +figure(3) +loglog(1:200,error4); +hold on +x = 1:200; +loglog(x,1./x) +hold on +loglog(x,1./sqrt(x)) +legend('Error of g', '1/N', '1/sqrt(N)') +title('Error for part 4') + + +%% +%Functions + +% Function to generate sinusoidal input signal +function u = generateSinusoidalInput(samples, periods, sinusoids) + omega = 2*pi/samples*round(samples/2*(1:sinusoids)'/(sinusoids+1)); % new + + thetas = pi*(1:sinusoids)'.*((1:sinusoids)'-1)/sinusoids; % Schroeder phase + + + k = (0:samples-1)'; + u = zeros(samples,1); + for i = 1:sinusoids + u = u + 2*cos(k*omega(i)+thetas(i)); + end + + u = 2*u /norm(u,"inf"); + + u = repmat(u,periods,1); + + +end + + +function y = generateOutput(u,noise,a,b,c, periods, samples, G) + y = zeros(periods*samples,1); + y_neg1 = 0; + u_neg1 = 0; + u_neg2 = 0; + y(1,1) = a*y_neg1 + b*u_neg1 + c*u_neg2 + noise(1,1); + y(2,1) = a*y(1) +b*u(1) +c*u_neg2 + noise(2,1); + y_neg1 = y(2); + u_neg1 = u(2); + u_neg2 = u(1); + for i = 3:length(u) + y(i) = a*y_neg1 +b*u_neg1 +c*u_neg2; + y_neg1 = y(i); + u_neg1 = u(i); + u_neg2 = u(i-1); + end + + % y = lsim(G,u); % Better this than the above one + + y = y + noise; + + + +end + + +function g = timeEstimation(u, tau_max, y) + c_toep = u; + r = [u(1), zeros(1,tau_max)]; + phi = toeplitz(c_toep,r); + + g = phi\y; + % other approach + % phy2 = zeros(length(u), tau_max); + % for i = 1:tau_max+1 + % phy2(i:end,i) = u(1:end-i+1); + % end + % g2 = phy2\y; + +end diff --git a/code/week08/SystemID_HS2023_Solution06.pdf b/code/week08/SystemID_HS2023_Solution06.pdf new file mode 100644 index 0000000..2c45bf2 Binary files /dev/null and b/code/week08/SystemID_HS2023_Solution06.pdf differ diff --git a/code/week09/Ex7_yhofmann.m b/code/week09/Ex7_yhofmann.m new file mode 100644 index 0000000..d84a6b0 --- /dev/null +++ b/code/week09/Ex7_yhofmann.m @@ -0,0 +1,136 @@ +%% +%1) +clear all +close all + +sigma_u = 1; +sigma_e = 0.1; +number_u = 200; +A = [1 0.25 -0.2 0.1 0.05]; +B = [0.6 0.3 -0.05 0]; +G = tf(B,A,1); +H = tf(1,A,1); + +u_id = sigma_u * randn(number_u,1); +e_id = sigma_e *randn(number_u,1); +v_id = lsim(H,e_id); +y_id = lsim(G,u_id) + v_id; + +%% +% 2.) +n_a = 1:8; +n_b = n_a; +system = cell(2,8); + + +% implement ARX +for i = 1:length(n_a) + + [system{1,i}, system{2,i}] = systemARX(i,i,y_id,u_id); + +end + + +%% +% 3.) +u_val = sigma_u * randn(number_u,1); +e_val = sigma_e * randn(number_u,1); + +y_val= lsim(G,u_val) + lsim(H,e_val); + +%% +% 4.) + +% Predict output for each and plot +y_pred_val = cell(1,8); +for i =1:length(n_a) + func = tf(system{2,i}',[1,system{1,i}'],1); + y_pred_val{i} = lsim(func,u_val); + figure(i) + plot(y_pred_val{i}); + hold on + plot(y_val); + title(sprintf('Predicted for n_a=n_b = %d',i)) + legend('Predicted', 'True'); +end + + +% calculate for identification +y_pred_id = cell(1,8); +for i =1:length(n_a) + func = tf(system{2,i}',[1,system{1,i}'],1); + y_pred_id{i} = lsim(func,u_id); +end + +%% +% 5.) + +% Compute error for each and plot +error_val = zeros(1,8); +error_id = zeros(1,8); +for i = 1:length(n_a) + error_val(1,i) = mean((y_val - y_pred_val{1,i}).^2); + error_id(1,i) = mean((y_id - y_pred_id{1,i}).^2); + +end + +figure(9) +plot(n_a,error_val); +hold on +plot(n_a, error_id); +xlabel('Order ARX') +title('Error') +legend('Validation', 'Identification') + + +%% +%6.) + +% Compute residuals for each and plot +residual = zeros(number_u,8); +for i = 1: length(n_a) + residual(:,i) = (y_val - y_pred_val{1,i}); +end +figure(100) +plot(1:number_u,residual); +title('Residual for different order') +legend('ARX1', 'ARX2', 'ARX3', 'ARX4', 'ARX5', 'ARX6', 'ARX7', 'ARX8') + + +%% +%7.) + +% Compare freq response +figure (200) +bode(G) +hold on + +for i = 1:length(n_b) + bode(tf(system{2,i}',[1, system{1,i}'],1),u_val) +end +hold off +legend('True', 'ARX1', 'ARX2', 'ARX3', 'ARX4', 'ARX5', 'ARX6', 'ARX7', 'ARX8') + +%% +%Functions +function [A,B] = systemARX(orderA, orderB, y, u) % can actually assume that was 0 before so no truncation for u and y needed +phi = zeros(length(y),orderA+orderB); +for i = 1:orderA + y_temp = -[zeros(i,1);y(1:end-i)]; + phi(:,i) = y_temp; +end + +for j = 1:orderB + u_temp =[zeros(j,1);u(1:end-j)]; + phi(:,orderA + j) = u_temp; +end + +theta = phi\y; +A = theta(1:orderA); +B = theta(orderA +1:end); + + + +%could also calculte G here and then make this the return value +%G = tf(B', [1, A'], 1); +end \ No newline at end of file diff --git a/code/week09/SystemID_HS2023_Solution07.pdf b/code/week09/SystemID_HS2023_Solution07.pdf new file mode 100644 index 0000000..df814bc Binary files /dev/null and b/code/week09/SystemID_HS2023_Solution07.pdf differ diff --git a/code/week10/Ex8_yhofmann.m b/code/week10/Ex8_yhofmann.m new file mode 100644 index 0000000..467e2c0 --- /dev/null +++ b/code/week10/Ex8_yhofmann.m @@ -0,0 +1,175 @@ +clear all, close all +%% +%a) +[u, y] = HS2023_SysID_Exercise_08_GenerateData(20914032); + +%% +%b) + +r = 20; +g_LSQ = timeEstimation(u,r,y); + + + +%% +%c) +alpha = 0.5; +gamma = 100; + +P = kernelTC(r,alpha); +phi = calculatePhi(u,r); +g_hat = (phi'*phi+ gamma*TCKernel_inv(alpha,r))\(phi'*y); + + + +%% +%d) +u_train = u(1:0.7*length(u)); +u_test = u(0.7*length(u)+1:end); + +y_train = y(1:0.7*length(y)); +y_test = y(0.7*length(y)+1:end); + + +%% +%e) +%TODO +r_candidates = 2:1:30; +alpha_candidates = 0:0.05:0.95; +find_best_alpha_r = zeros(length(alpha_candidates),length(r_candidates),2); +for i = 1:length(alpha_candidates) + find_best_alpha_r(i,:,1) = r_candidates; +end +for j = 1:length(r_candidates) + find_best_alpha_r(:,j,2) = alpha_candidates'; +end + + +tuned_alpha = 0; +tuned_r = 2; +tuned_g_hat = cell(1,1); +smallest_error = 10e20; + +error_ = zeros(length(alpha_candidates),length(r_candidates)); +%g_hat_tuned = cell(length(alpha_candidates),length(r_candidates)); +%y_hat_tuned = zeros(length(u_train),size(alpha_candidates) ) + +for i= 1:length(alpha_candidates) + for j = 1:length(r_candidates) + alpha = find_best_alpha_r(i,j,2); + r = find_best_alpha_r(i,j,1); + + P = kernelTC(r,alpha); + phi = calculatePhi(u_train,r); + % g_hat_tuned{i,j} = (phi'*phi+ gamma*inv(P))\(phi'*y); + g_hat_current = (phi'*phi + gamma*TCKernel_inv(alpha,r))\(phi'*y_train); + y_estimated = pulseResponseOutput(u_test,g_hat_current); + error_(i,j) = mean((y_estimated-y_test).^2); % debug + if(mean((y_estimated-y_test).^2) < smallest_error) + smallest_error = mean((y_estimated-y_test).^2) + tuned_r = r; + tuned_alpha = alpha; + tuned_g_hat{1,1} = g_hat_current; + + if(smallest_error) < 1 + figure(21436) + plot(y_estimated) + hold on + plot(y_test) + legend('estimated', 'true') + i = 2 + close all + end + end + end +end + +disp("Best alpha = " + num2str(tuned_alpha)); +disp("Best r = " + num2str(tuned_r)); + + + + + +%% +%f) +%TODO + +figure(100) +stem(g_LSQ) +hold on +stem(g_hat) +hold on +stem(tuned_g_hat{1,1}) + +title('Different pulse responses') +xlabel('parameter number') +ylabel('Value') +legend('Least square','Regularized proposed', 'Regularized tuned') + + + +figure(200) +surf(error_) +xlabel('order r') +ylabel('alpha') +zlabel('error') +title('Error for different parameter values') + +%% FUNCTIONS +function g = timeEstimation(u, tau_max, y) + c_toep = u; + r = [u(1), zeros(1,tau_max)]; + phi = toeplitz(c_toep,r); + + + + g = phi\y; + % other approach + % phy2 = zeros(length(u), tau_max); + % for i = 1:tau_max+1 + % phy2(i:end,i) = u(1:end-i+1); + % end + % g2 = phy2\y; + +end + + +function P = kernelTC(N,alpha) +P = zeros(N,N); +for i = 1:N + for j = 1:N + P(i,j) = alpha^max(i,j); + end +end +end + + +function Phi = calculatePhi(u,N) + c_toep = u; + r = [u(1), zeros(1,N-1)]; + Phi = toeplitz(c_toep,r); + + + +end + + + +function y = pulseResponseOutput(u,pulse) + tau_max = length(pulse) -1; + c_toep = u; + r = [u(1), zeros(1,tau_max)]; + phi = toeplitz(c_toep,r); + + y = phi * pulse; + +end + + +function invP = TCKernel_inv(alpha, n) +diag0 = alpha .^ (-(1:n)') / (1-alpha); +diag0(2:n-1) = diag0(2:n-1) * (1+alpha); +diag1 = -alpha .^ (-(1:n-1)') / (1-alpha); +invP = diag(diag0, 0) + diag(diag1, +1) + diag(diag1, -1); +end \ No newline at end of file diff --git a/code/week10/SystemID_HS2023_Exercise8__regularization_solution.pdf b/code/week10/SystemID_HS2023_Exercise8__regularization_solution.pdf new file mode 100644 index 0000000..1699851 Binary files /dev/null and b/code/week10/SystemID_HS2023_Exercise8__regularization_solution.pdf differ diff --git a/code/week11/Ex9_yhofmann.m b/code/week11/Ex9_yhofmann.m new file mode 100644 index 0000000..8f73381 --- /dev/null +++ b/code/week11/Ex9_yhofmann.m @@ -0,0 +1,120 @@ +clear all, close all + +N = 10^3; +A = [1 -1.5 0.7]; +B = [0 1 0.4]; +C = [1 -1.1 0.4]; + + +e = randn(N,1); +u = getU(N,e); +y_true = lsim(tf(B,A,1),u) + lsim(tf(C,A,1),e); +y_true_filtered = lsim(tf(1,C,1),y_true); +u_filtered = lsim(tf(1,C,1),u); + +%% +%1.) +[A_LS, B_LS] = LS_from_ARMAX(N,y_true_filtered,u_filtered); + +phi = zeros(N,4); +phi(2,:) = [-y_true_filtered(1) 0 u_filtered(1) 0]; +for k = 3:N + phi(k,:) = [-y_true_filtered(k-1) -y_true_filtered(k-2) u_filtered(k-1) u_filtered(k-2)]; +end + +theta = phi\y_true_filtered; + + +%% +%2.) +e_val = randn(N,1); +u_val = getU(N,e_val); +%u_val_filtered = lsim(1/C,u_val); + +y_val = lsim(tf(B,A,1),u_val) + lsim(tf(C,A,1),e_val); +y_pred = lsim(tf(B_LS,A_LS,1), u_val) + lsim(tf(C,A_LS,1),e_val); + +figure(1) +plot(y_val) +hold on +plot(y_pred) +legend('True', 'Predicted') +xlabel('step k') +ylabel('value') + + + +%% +%3.) + +num_Exp = 1000; +A_3 = zeros(3,num_Exp); +B_3 = zeros(3,num_Exp); + + +for k = 1:num_Exp + e_exp = randn(N,1); + u_exp = getU(N,e_exp); + + y_exp = lsim(tf(B,A,1),u_exp) + lsim(tf(C,A,1),e_exp); + y_exp_filtered = lsim(tf(1,C,1),y_exp); + u_exp_filtered = lsim(tf(1,C,1),u_exp); + + [A_3(:,k), B_3(:,k)] = LS_from_ARMAX(N,y_exp_filtered,u_exp_filtered); + + +end + + +figure(2) +sgtitle('Histograms LS estimates different noise realizations (n = 1000)') + +subplot(2,2,1) +histogram(A_3(2,:)) +title('a_1') + +subplot(2,2,2) +histogram(A_3(3,:)) +title('a_2') + +subplot(2,2,3) +histogram(B_3(2,:)) +title('b_1') + + +subplot(2,2,4) +histogram(B_3(3,:)) +title('b_2') + + + + + + +%% +%FUNCTIONS + +function u = getU(N, e) + u = zeros(N,1); + u(2) = 0.9*e(1); + for k = 3:N + u(k) = 0.14 *u(k-1) +0.12*u(k-2) + 0.9 * e(k-1) + 0.3 * e(k-2); + end +end + + +function [A,B] = LS_from_ARMAX(N,y,u) + +phi = zeros(N,4); +phi(2,:) = [-y(1) 0 u(1) 0]; +for k = 3:N + phi(k,:) = [-y(k-1) -y(k-2) u(k-1) u(k-2)]; +end + +theta = phi\y; + +B = [0, theta(3) , theta(4)]; +A = [1 , theta(1) , theta(2)]; + +end + diff --git a/code/week11/SystemID_HS2023_Solution09.pdf b/code/week11/SystemID_HS2023_Solution09.pdf new file mode 100644 index 0000000..2e4c2c2 Binary files /dev/null and b/code/week11/SystemID_HS2023_Solution09.pdf differ diff --git a/code/week12/Ex10_yhofmann.m b/code/week12/Ex10_yhofmann.m new file mode 100644 index 0000000..4c70401 --- /dev/null +++ b/code/week12/Ex10_yhofmann.m @@ -0,0 +1,60 @@ +clear all, close all +%% Load and initialize +data = importdata("p3_data.mat"); +r1 = data.p3_r1; +y1 = data.p3_y1; + +z = tf('z'); +F = 0.5*z^-2 + z^-1; +F2 = tf([0, 1,0.5],1,-1, 'Variable','z^-1'); % also like this to avoid a fractionn +%% +%1) +u1 = r1 - lsim(F,y1); + +[A,B] = systemARX_no_Initial_Conditions(2,1,y1,u1); +a1 = A(1) +a2 = A(2) +b1 = B + +G_est = tf([b1 0],[1 a1 a2],1); +Est_closed_system = feedback(G_est,F); +y_est = lsim(Est_closed_system,r1); + +figure(1) +plot(y1); +hold on +plot(y_est) +title('Direct Method') +legend('Data', 'Estimated') +grid on + + + +%% +% 2.) +% MATLAB code would be exactly the same as above but error term (which does +% not matter for calculation is different and is now be unbiased) + + +%% +%Functions + +%!!!!!! Instead of only writng it for the truncated case can also use the +%normal one and then simply for the \ solving disregard the first 3 rows +function [A,B] = systemARX_no_Initial_Conditions(orderA, orderB, y, u) % here cannot assume that initial conditition was 0 +phi = zeros(length(y) - orderA -orderB+1,orderA+orderB); +for i = 1:orderA + y_temp = -[y(orderA+orderB - i:end-i)]; + phi(:,i) = y_temp; +end + +for j = 1:orderB + u_temp =[u(orderA + orderB -j:end-j)]; + phi(:,orderA + j) = u_temp; +end + +theta = phi\y(orderA+orderB:end); +A = theta(1:orderA); +B = theta(orderA +1:end); + +end \ No newline at end of file