# offline ( Gauss Newtone approach ) Prediction error method for ARMAX Model, Algorithm Implemented but not good results

3 ビュー (過去 30 日間)
Reda Chtiba 2022 年 5 月 12 日

% 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
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)
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
theta=theta+0.8*P;%GaussNewtone Iterations
end

サインインしてコメントする。

### 回答 (1 件)

Aman 2024 年 2 月 7 日
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!

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeOil, Gas & Petrochemical についてさらに検索

R2022a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by