necesito ayuda con esto
Este es un programa que he hecho para identificacion parametrica del sistema. El metodo tiene que ser NEwton Raphson con funciones de sensibilidad, pero la solucion diverge y no logro solucionarlo. Me podeis ayudar?
clear all
close all
%visualisation des jeux de données
xm=[1 1 1 1 1 1 1 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 -0.25 -0.25 -0.25 -0.25 -0.25 -0.25 -0.25 -0.25 -0.25 0.5 0.5 0.5 0.5 0.5 0.5];
ym=[0.043 0.13 0.23 0.36 0.46 0.55 0.62 0.66 0.66 0.68 0.6 0.59 0.545 0.59 0.55 0.53 0.40 0.34 0.24 0.12 0.14 0.04 -0.01 -0.04 -0.06 0.00 0.10 0.17 0.21 0.26];
subplot(2,1,1);
plot(xm,'b'),title('entree mesuree')
subplot(2,1,2);
plot(ym,'r'),title('sortie mesuree');pause
%initialisation des paramètres
n=length(xm);
a=0.8;
b=0.2;
%trace de la reponse initiale
y=zeros(1,n);
for j=1:n-1
y(j+1)=a*y(j)+b*xm(j);
end
figure(2)
j=(1:30);
plot(j,ym,'o',j,y),title('Evolution de la reponse identifiee');pause
%Calcul de l'erreur, du gradient, des f. de sensibilite et du hessian
p(:,1)=[a;b];
sa=zeros(1,n);
sb=zeros(1,n);
L=0.8;
ite=2;
%do{
for ite=2:10
G=zeros(2,1);
H=zeros(2,2);
S=[1 0;0 1];
e=y-ym;
for k=1:n-1
%y(k+1)=a*y(k)+b*xm(k); %calcul de la sortie du modele
sa(k+1)=a*sa(k)+y(k); %equations de sensibilite
sb(k+1)=a*sb(k)+xm(k);
end
for k=1:n-1
G(1,1)=G(1,1)+e(k)*sa(k); %components du gradient
G(2,1)=G(2,1)+e(k)*sb(k);
H(1,1)=H(1,1)+(sa(k)^2)*L^(n-k); %components du hessian
H(1,2)=H(1,2)+sa(k)*sb(k)*L^(n-k);
H(2,2)=H(2,2)+(sb(k)^2)*L^(n-k);
end
H(2,1)=H(1,2);
H=H+1e-3*S;
p(:,ite)=p(:,ite-1)-(L^2)*inv(H)*G; %actualisation des parametres
a=p(1,ite);
b=p(2,ite);
%trace de la reponse
i=(1:n);
for j=1:n-1
y(j+1)=a*y(j)+b*xm(j);
end
figure(3)
plot(i,ym,'x',i,y),title('evolution de la reponse');pause
end
%}while norm(G)>1e-4
%trace de resultats
j=(1:ite);
figure(4)
subplot(3,1,1)
plot(i,ym,'x',i,y),title('sortie du modele');
subplot(3,1,2)
plot(j,p(1,:)),title('parametre a');
subplot(3,1,3)
plot(j,p(2,:)),title('parametre b');
clear all
close all
%visualisation des jeux de données
xm=[1 1 1 1 1 1 1 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 -0.25 -0.25 -0.25 -0.25 -0.25 -0.25 -0.25 -0.25 -0.25 0.5 0.5 0.5 0.5 0.5 0.5];
ym=[0.043 0.13 0.23 0.36 0.46 0.55 0.62 0.66 0.66 0.68 0.6 0.59 0.545 0.59 0.55 0.53 0.40 0.34 0.24 0.12 0.14 0.04 -0.01 -0.04 -0.06 0.00 0.10 0.17 0.21 0.26];
subplot(2,1,1);
plot(xm,'b'),title('entree mesuree')
subplot(2,1,2);
plot(ym,'r'),title('sortie mesuree');pause
%initialisation des paramètres
n=length(xm);
a=0.8;
b=0.2;
%trace de la reponse initiale
y=zeros(1,n);
for j=1:n-1
y(j+1)=a*y(j)+b*xm(j);
end
figure(2)
j=(1:30);
plot(j,ym,'o',j,y),title('Evolution de la reponse identifiee');pause
%Calcul de l'erreur, du gradient, des f. de sensibilite et du hessian
p(:,1)=[a;b];
sa=zeros(1,n);
sb=zeros(1,n);
L=0.8;
ite=2;
%do{
for ite=2:10
G=zeros(2,1);
H=zeros(2,2);
S=[1 0;0 1];
e=y-ym;
for k=1:n-1
%y(k+1)=a*y(k)+b*xm(k); %calcul de la sortie du modele
sa(k+1)=a*sa(k)+y(k); %equations de sensibilite
sb(k+1)=a*sb(k)+xm(k);
end
for k=1:n-1
G(1,1)=G(1,1)+e(k)*sa(k); %components du gradient
G(2,1)=G(2,1)+e(k)*sb(k);
H(1,1)=H(1,1)+(sa(k)^2)*L^(n-k); %components du hessian
H(1,2)=H(1,2)+sa(k)*sb(k)*L^(n-k);
H(2,2)=H(2,2)+(sb(k)^2)*L^(n-k);
end
H(2,1)=H(1,2);
H=H+1e-3*S;
p(:,ite)=p(:,ite-1)-(L^2)*inv(H)*G; %actualisation des parametres
a=p(1,ite);
b=p(2,ite);
%trace de la reponse
i=(1:n);
for j=1:n-1
y(j+1)=a*y(j)+b*xm(j);
end
figure(3)
plot(i,ym,'x',i,y),title('evolution de la reponse');pause
end
%}while norm(G)>1e-4
%trace de resultats
j=(1:ite);
figure(4)
subplot(3,1,1)
plot(i,ym,'x',i,y),title('sortie du modele');
subplot(3,1,2)
plot(j,p(1,:)),title('parametre a');
subplot(3,1,3)
plot(j,p(2,:)),title('parametre b');