commit
2a25844091
10 changed files with 414 additions and 0 deletions
@ -0,0 +1,5 @@ |
|||
load('./SysID_Exercise_1.mat') |
|||
|
|||
w = [x.' ; 2.*x.'.^2-1].'; |
|||
|
|||
theta = inv(w.' * w) * w.'*(y-0.2) |
|||
@ -0,0 +1,11 @@ |
|||
% Ridge regression |
|||
|
|||
load('./SysID_Exercise_1.mat') |
|||
|
|||
w = [x.' ; 2.*x.'.^2-1].'; |
|||
|
|||
u = [1; 0.4]; |
|||
sigma_p = 0.01; |
|||
sigma = 0.1; |
|||
|
|||
theta = inv(1/sigma * w.' * w + 1/sigma_p*eye(2)) * (1/sigma*w.'*(y-0.2)+u/sigma_p) |
|||
@ -0,0 +1,17 @@ |
|||
load('./SysID_Exercise_1.mat') |
|||
|
|||
w = [x.' ; 2.*x.'.^2-1].'; |
|||
|
|||
theta_ml = inv(w.' * w) * w.'*(y-0.2); |
|||
|
|||
u = [1; 0.4]; |
|||
sigma_p = 0.01; |
|||
sigma = 0.1; |
|||
|
|||
theta_map = inv(1/sigma * w.' * w + 1/sigma_p*eye(2)) * (1/sigma*w.'*(y-0.2)+u/sigma_p); |
|||
|
|||
y_hat_ml = theta_ml(1) * x_v + theta_ml(2) * (2.*x_v.^2-1)+0.2; |
|||
y_hat_map = theta_map(1) * x_v + theta_map(2) * (2.*x_v.^2-1)+0.2; |
|||
|
|||
MSE_MLE = sum((y_hat_ml - y_v).^2) / length(y_v) |
|||
MSE_MAP = sum((y_hat_map - y_v).^2) / length(y_v) |
|||
@ -0,0 +1,5 @@ |
|||
function [ret] = eccentricity(a, b, c, d, e, f) |
|||
matrix = [a, b/2, d/2; b/2, c, e/2; d/2, e/2,f]; |
|||
eta = - sign(det(matrix)); |
|||
ret = sqrt((2*sqrt((a-c)^2+b^2)) / (eta*(a+c)+sqrt((a-c)^2+b^2))); |
|||
end |
|||
@ -0,0 +1,20 @@ |
|||
% Coment out the comet not to compute |
|||
load('./C-2017-K2.mat') % First comet |
|||
% load('./P-2011-NO1.mat') % Second comet |
|||
|
|||
x = xy(:,1); |
|||
y = xy(:,2); |
|||
|
|||
X = [x.^2, x.*y, y.^2, x, y]; |
|||
Y = -ones(size(X, 1), 1); |
|||
|
|||
theta = X.' * X \ X.' * Y |
|||
|
|||
eccentricity(theta(1), theta(2), theta(3), theta(4), theta(5), 1) |
|||
|
|||
syms xp yp |
|||
fimplicit(theta(1)*xp^2 + theta(2)*xp*yp + theta(3)*yp^2 + theta(4)*xp + theta(5)*yp == 1) |
|||
hold on |
|||
|
|||
scatter(x, y, 'r'); |
|||
|
|||
@ -0,0 +1,6 @@ |
|||
% Assume that u is periodic |
|||
N=400; |
|||
t = 0:N-1; |
|||
u = 1.*randn(N, 1); |
|||
|
|||
lags = [0:(N-1)/2]; |
|||
@ -0,0 +1,33 @@ |
|||
%% 1. Generate e(k), random noise |
|||
N = 1024; |
|||
e = randn(N, 1); |
|||
|
|||
%% 2. Calculate and plot the periodogram of e |
|||
% Fourier transform |
|||
E = fft(e); |
|||
omega = (2*pi/N)*[0:N-1]'; |
|||
% Positive frequencies |
|||
pos_idx = find(omega >= 0.0 & omega <= pi); |
|||
% Plot |
|||
E_periodogram = abs(E(pos_idx)).^2/N; |
|||
loglog(omega(pos_idx), E_periodogram) |
|||
|
|||
%% 3. Given plant P, calculate the response (w) |
|||
num = [1 0.5]; |
|||
den = [1 0 0.25]; |
|||
P = tf(num, den, 1.0); |
|||
w = lsim(P, e); |
|||
|
|||
%% 4. Calculate the periodogram of w(k) |
|||
W = fft(w); |
|||
W_periodogram = abs(W(pos_idx)).^2/N; |
|||
loglog(omega(pos_idx), W_periodogram) |
|||
hold on |
|||
|
|||
%% 5. Compare periodogram w to square of Bode Plot |
|||
[mag, phase, wout] = bode(P, omega(pos_idx)); |
|||
mag = squeeze(mag); |
|||
loglog(omega(pos_idx), abs(mag).^2) |
|||
|
|||
diff = abs(W_periodogram - mag); |
|||
loglog(omega(pos_idx), diff) |
|||
@ -0,0 +1,59 @@ |
|||
G = tf([0.0565 0.1013], [1 -1.4183 1.5894 -1.3161 0.8864], 1.0); |
|||
|
|||
%% a) Generate sets of identification data |
|||
N = 200; |
|||
|
|||
u_idin = idinput(N); |
|||
u_unif = 2 * rand(N, 1) - 1; |
|||
|
|||
% for 3c) |
|||
N = 5 * N; |
|||
u_idin = repmat(u_idin, 5, 1); |
|||
u_unif = repmat(u_unif, 5, 1); |
|||
|
|||
e = 0.1 * randn(N, 1); |
|||
y_idin = lsim(G, u_idin) + e; |
|||
y_unif = lsim(G, u_unif) + e; |
|||
|
|||
%% b) Estimate G(z) by EFTE |
|||
U_idin = abs(fft(u_idin)); |
|||
U_unif = abs(fft(u_unif)); |
|||
|
|||
Y_idin = abs(fft(y_idin)); |
|||
Y_unif = abs(fft(y_unif)); |
|||
|
|||
G_idin = Y_idin ./ U_idin; |
|||
G_unif = Y_unif ./ U_unif; |
|||
|
|||
omega = (2*pi/N)*[0:N-1]; |
|||
pos_idx = find(omega >= 0.0 & omega <= pi); |
|||
|
|||
H = abs(squeeze(freqresp(G, omega(pos_idx)))); |
|||
|
|||
err_idin = abs(H - G_idin(pos_idx)); |
|||
err_unif = abs(H - G_unif(pos_idx)); |
|||
|
|||
subplot(2, 1, 1) |
|||
loglog(omega(pos_idx), H, 'green'); |
|||
hold on |
|||
loglog(omega(pos_idx), G_idin(pos_idx), '--b'); |
|||
loglog(omega(pos_idx), G_unif(pos_idx), '-.r'); |
|||
hold off |
|||
legend('True plant', 'EFTE random binary', 'EFTE uniform distribution') |
|||
|
|||
subplot(2, 1, 2) |
|||
loglog(omega(pos_idx), err_idin, '--b') |
|||
hold on |
|||
loglog(omega(pos_idx), err_unif, '-.r') |
|||
hold off |
|||
legend('EFTE random binary error', 'EFTE uniform distribution error') |
|||
|
|||
fprintf('err_idin total: %f\n', mean(err_idin)) |
|||
fprintf('err_unif total: %f\n', mean(err_unif)) |
|||
|
|||
%% c) |
|||
% 1) Obvious. Slightly better results |
|||
% 2) |
|||
% 3) Best solution. Discard first period. Average the rest and apply fft after |
|||
% y5 = mean(reshape(yp(N+1:end), N, []), 2) |
|||
% Y5_N = fft(y5) |
|||
@ -0,0 +1,163 @@ |
|||
function [R_u,lags,omegas,U_cos,U_cos_Wf,U_cos_Wt] = HS2022_SysID_Exercise_05_12345678() |
|||
|
|||
%% Output format specification |
|||
% R_u must be a 3xN matrix |
|||
% lags must be a 1xN vector |
|||
% omegas must be a 1xN vector |
|||
% U_cos must be a 1xN vector |
|||
% U_cos_Wf must be a 1xN vector |
|||
% U_cos_Wt must be a 1xN vector |
|||
%% Generate data |
|||
|
|||
% Extract Legi from Filename |
|||
name=mfilename; |
|||
LegiNumber= name(end-7:end); |
|||
|
|||
[u_prbs,u_rand,u_cos] = HS2022_SysID_Exercise_05_GenerateData(LegiNumber); |
|||
|
|||
N = length(u_prbs); |
|||
|
|||
%% General instructions for solution |
|||
|
|||
% Change the filename of this function, both in the function definition |
|||
% above and in the filename in the folder |
|||
|
|||
% Use the input signals u_prbs, u_randn and u_cos to solve the problem. |
|||
|
|||
% Modify your code in the next sections, and return the variables |
|||
% requested. |
|||
|
|||
% If you skip one part of the problem, return the empty vectors as already |
|||
% provided in the code |
|||
|
|||
%% 1. Calculation of autocorrelation |
|||
|
|||
R_u = NaN * ones(3,N); |
|||
|
|||
lags = NaN * ones(1,N); |
|||
|
|||
|
|||
%% 2. Smoothing |
|||
|
|||
omegas = NaN * ones(1,N); |
|||
|
|||
U_cos = NaN * ones(1,N); |
|||
|
|||
U_cos_Wf = NaN * ones(1,N); |
|||
|
|||
U_cos_Wt = NaN * ones(1,N); |
|||
|
|||
|
|||
%% Functions |
|||
|
|||
function [lags_w,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_w = [floor(-N/2+1):floor(N/2)]'; |
|||
wHann = 0*lags_w; |
|||
idx = find(abs(lags_w) <= gamma); |
|||
wHann(idx) = 0.5*(1+cos(pi*lags_w(idx)/gamma)); |
|||
|
|||
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_w = [floor(-N/2+1):floor(N/2)]'; |
|||
wHann = 0*lags_w; |
|||
idx = find(abs(lags_w) <= gamma); |
|||
wHann(idx) = 0.5*(1+cos(pi*lags_w(idx)/gamma)); |
|||
|
|||
% |
|||
zidx = find(lags_w==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_w; |
|||
|
|||
end |
|||
|
|||
end |
|||
@ -0,0 +1,95 @@ |
|||
a=5/8; |
|||
b=11/8; |
|||
cv=0.05; |
|||
|
|||
%% Exercise 6.1 |
|||
N_p = 2; |
|||
N = 2; |
|||
tau_max = 80; |
|||
u_max = 2; |
|||
|
|||
g61 = estimate_g(N, N_p, tau_max, u_max, a, b, cv, true, true); |
|||
|
|||
%% Exercise 6.2 |
|||
|
|||
N_p = 1; |
|||
N = 1; |
|||
tau_max = 80; |
|||
u_max = 2; |
|||
g0 = estimate_g(N, N_p, tau_max, u_max, a, b, cv, false, true); |
|||
|
|||
err = zeros(30, 1); |
|||
for N_p=2:30 |
|||
g = estimate_g(N, N_p, tau_max, u_max, a, b, cv, true, true); |
|||
err(N_p) = norm(g - g0); |
|||
end |
|||
|
|||
plot(2:30, err(2:end), 'b'); |
|||
xlabel("Number of Periods") |
|||
ylabel("||g_i - g_0||_2") |
|||
grid on |
|||
hold on |
|||
|
|||
%% Exercise 6.3 |
|||
|
|||
N = 1; |
|||
tau_max = 80; |
|||
u_max = 2; |
|||
err = zeros(30, 1); |
|||
for N_p=2:30 |
|||
g = estimate_g(N, N_p, tau_max, u_max, a, b, cv, true, false); |
|||
err(N_p) = norm(g - g0); |
|||
end |
|||
|
|||
plot(2:30, err(2:end)) |
|||
|
|||
%% Exercise 6.4 |
|||
|
|||
N_p = 2; |
|||
tau_max = 80; |
|||
u_max = 2; |
|||
|
|||
err = zeros(30, 1); |
|||
for N=2:30 |
|||
g = estimate_g(N, N_p, tau_max, u_max, a, b, cv, true, true); |
|||
err(N) = norm(g - g0); |
|||
end |
|||
|
|||
plot(2:30, err(2:end)) |
|||
|
|||
|
|||
%% Exercise function |
|||
function g=estimate_g(N, N_p, tau_max, u_max, a, b, cv, noise, PRBS) |
|||
g = zeros(tau_max, 1); |
|||
for j=1:N |
|||
% Compute the noise vector |
|||
if noise |
|||
v = sqrt(cv) * randn(N_p * tau_max, 1); |
|||
else |
|||
v = zeros(N_p * tau_max, 1); |
|||
end |
|||
|
|||
% Compute the input vector |
|||
if PRBS |
|||
u = u_max * idinput(tau_max, 'PRBS'); |
|||
else |
|||
u = u_max * (rand(tau_max, 1) - 0.5) * 2.0; |
|||
end |
|||
u = repmat(u, N_p, 1); |
|||
|
|||
% Compute the outputs of the system |
|||
y = zeros(size(u)); |
|||
y(1) = v(1); |
|||
for i=2:size(y, 1) |
|||
y(i) = a * y(i-1) + b * u(i-1) + v(i); |
|||
end |
|||
|
|||
% Build the Phi matrix |
|||
% phi = [toeplitz(y, zeros(tau_max, 1)), toeplitz(u, zeros(tau_max, 1))]; |
|||
phi = toeplitz(u, zeros(tau_max, 1)); |
|||
|
|||
% Solve for g |
|||
g = g + phi \ y; |
|||
end |
|||
g = g / N; |
|||
end |
|||
Loading…
Reference in new issue