Now suppose you are asked to develop an online wind-speed prediction algorithm t

Question

Now suppose you are asked to develop an online wind-speed prediction algorithm that predicts wind
speed 10 minutes into the future. Suppose that you decide to use an AR model of the following form:

Part A

Data collection and preparation:
Go to http://wind.nrel.gov/Web_nrel/ , enter “1” into the form, and download the wind data of 2006. The data is stored in a “.csv” file. Suppose that you name this file as “wind.csv”. In matlab, you can read in the data using ”M = csvread(’wind.csv’,2,1)”. The first column of matrix M will be the wind speed data (i.e., the measurement data y(1), y(2), . . .) we need for this problem.

Part B

I am really looking for part b

So far all I have gotten is

M = csvread('wind.csv',2,1);
y = M(:,1);
y = y(1:3000);
n=(1:20);

I tried looking up ar and arx in matlab but I am stuck

Explanation / Answer

Solution :

%% Load data

datafile = 'pmpat010.txt';

data = load(datafile);

t=data(1,:); %[s]

Ca=data(2,:); %[mM] arterial input function

Ct=data(3,:); %[mM] tissue response

figure; plot(t,Ca,'.b',t,Ct,'.r'); hold on

xlabel('Time [s]'); ylabel('[mM]')

legend('AIF','Ct')

N=length(Ct);

%% Simulate model with initial parameter values

%Simulate

Ktrans = 2e-3; ve = 0.1;

p = [Ktrans, ve];

[y,t] = compart_lsim(t,Ca,p);

plot(t,y,'k');

%% Estimate parameters and variances

optim_options = optimset('Display', 'iter',...

...%'TolFun', 1e-6,... %default: 1e-4

...%'TolX', 1e-6,... %default: 1e-4

'LevenbergMarquardt', 'on'); %default: 'off'

%optim_options = [];

p0 = [Ktrans, ve];

[p,resnorm,residual,exitflag,OUTPUT,LAMBDA,Jacobian] = lsqnonlin(@compart_error, p0, [],[],optim_options, t,Ca,Ct);

disp(' ')

p

%Estimated parameter variance

Jacobian = full(Jacobian); %lsqnonlin returns the Jacobian as a sparse matrix

varp = resnorm*inv(Jacobian'*Jacobian)/N

stdp = sqrt(diag(varp)); %standard deviation is square root of variance

stdp = 100*stdp'./p; %[%]

disp([' Ktrans: ', num2str(p(1)), ' +/- ', num2str(stdp(1)), '%'])

disp([' ve: ', num2str(p(2)), ' +/- ', num2str(stdp(2)), '%']);

%% Simulate estimated model

[y,t] = compart_lsim(t,Ca,p);

figure; subplot(211); plot(t,Ct,'.r',t,y,'b');

xlabel('Time [s]'); ylabel('[mM]')

legend('data','model')

xi = Ct(:)-y(:); %same as residual from lsqnonlin

subplot(212); plot(t,xi)

xlabel('Time [s]'); legend('residuals xi')

assen=axis;

%% function to simulate compartment model

function [y,t] = compart_lsim(t,Ca,p)

Ktrans = p(1); ve = p(2);

num = [Ktrans*ve];

den = [ve, Ktrans];

sys = tf(num,den);

[y,t] = lsim(sys,Ca,t);

end

%% function to calculate MLE error

function xi = compart_error(p, t,Ca,Ct)

[y,t] = compart_lsim(t,Ca,p);

xi = Ct(:)-y(:); %make sure both vectors are columns

%figure(1); plot(t,Ct,'.',t,y,'g'); hold off; drawnow %uncomment this line to show

%datafir for each optimization iteration (slows down the execution)

end

end %main function

Related Questions

Navigate

• Browse (All)
• Browse N
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.