offline ( Gauss Newtone approach ) Prediction error method for ARMAX Model, Algorithm Implemented but not good results
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
% Simulation of an ARMAX model y(t)-0.7y(t-1)=0.5u(t-1)+e(t)+0.6e(t)
np=251;
a0=[1 -0.7 ]; b0=[0 0.5 ]; c0=[1 0.6]; % Set parameter values
u=randn(np,1)*sqrt(1); % Generate white noise with variance =1
e=randn(np,1)*sqrt(1); % Generate white noise with variance = 1
y=filter(b0,a0,u)+filter(c0,a0,e); % generate output data : y=(B/A)*u+(C/A)*e
z=[y u]; % store our inpute output data
x=rpem(z,[1 1 1 0 0 1],'ff',1); % build-in predction error method
x(end,:);
%% %%
theta=[-0.6 ;0.4 ; 0.6];% initiliazing
grad=[]; % dummy
iter_max=10;
for iter=1:iter_max
theta1=theta(1);
theta2=theta(2);
theta3=theta(3);
pe=filter([1 theta1],[1 theta3],y)-filter([0 theta2],[1 theta3],u);% Predictive errors, for the armax is : pe=(A/C)y(t)-(B/C)u(t)
for i=2:251 % Gradient for armax is given by grad(t)=(1/C)(-y(t-1),u(t-1),pe(t-1)) .
grad(1,:)=[0 0 0];% dummy value at time 1, to be deleted later
grad(i,:)=filter(1,[1 theta3],[-y(i-1) u(i-1) pe(i-1)]);% we started at i=2 so we can calculate using past value for other data
grad(1,:)=[];% delete first dummy row , we get 250 data
end
pe(1,:)=[];% delete first row , to get 250 data so we can compare it with grad
P=(grad'*grad)\(grad'*pe); % Gauss-Newtone Formula
theta=theta+0.8*P;%GaussNewtone Iterations
end
0 commentaires
Réponses (1)
Aman
le 7 Fév 2024
Hi Reda,
As per my understanding, you are using the Gauss-Newton method of optimizing the prediction error, but the result is not converging.
The matrix that has been created is near singular, so in that case the delta is not getting computed properly. In this case, you need to add a damping factor so that the proper delta is calculated. Refer to the below code, where I have added a damping factor of "1e-3".
% Simulation of an ARMAX model y(t)-0.7y(t-1)=0.5u(t-1)+e(t)+0.6e(t)
np=251;
a0=[1 -0.7 ];
b0=[0 0.5 ];
c0=[1 0.6]; % Set parameter values
u=randn(np,1)*sqrt(1); % Generate white noise with variance = 1
e=randn(np,1)*sqrt(1); % Generate white noise with variance = 1
y=filter(b0,a0,u)+filter(c0,a0,e); % generate output data : y=(B/A)*u+(C/A)*e
z=[y u]; % store our input output data
%%
y_prev = [0; y(1:end-1)]; % y(t-1)
u_prev = [0; u(1:end-1)]; % u(t-1)
e_prev = zeros(np, 1);
theta=[-0.6 ;0.4 ; 0.6];% initiliazing
iter_max=10;
error = [];
for iter=1:iter_max
residual = y - (theta(1)*y_prev + theta(2)*u_prev + theta(3)*e_prev);
J = [-y_prev, u_prev, e_prev];
lambda = 1e-3; % Damping factor
delta = ((J'*J + lambda*eye(size(J,2))) \ J') * residual;
theta = theta + delta;
e_prev = [0; residual(1:end-1)];
end
I hope it helps!
0 commentaires
Voir également
Catégories
En savoir plus sur Oil, Gas & Petrochemical dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!