# How do I get this code to converge?

3 visualizaciones (últimos 30 días)
Kaydian Quintero el 20 de Jul. de 2021
Respondida: Vaibhav el 12 de Abr. de 2024
% This code runs, but theta_hat should converge to 1, and e which is the
% error should converge to 0.
clc;
clear all;
close all;
T=120;
%Simulation step size
dt = 0.0001;
t(1)=0;
N=T/dt;
% Initial conditions
psi = [1/sqrt(2) 1/sqrt(2)]';
rho(:,:,1) = psi*psi';
psi_hat = [1/sqrt(5) 1/sqrt(5)]';
rho_hat(:,:,1) = psi_hat*psi_hat';
y(1) = 0;
y_hat(1) = 0;
theta_hat(1) = 1.5;
% Parameters definitions
A = 1;
w = 1;
omega=1;
theta = 1;
gamma = 0.1;
GAMMA = 1;
mu = [0 1;1 0];
P = [1 0;0 0];
H = [w/2 0;0 -w/2];
rho11(1)=rho(1,1,1);
rho12(1)=rho(1,2,1);
rho21(1)=rho(2,1,1);
rho22(1)=rho(2,2,1);
rho_hat11(1)=rho_hat(1,1,1);
rho_hat12(1)=rho_hat(1,2,1);
rho_hat21(1)=rho_hat(2,1,1);
rho_hat22(1)=rho_hat(2,2,1);
for k=1:N
t(k+1)=t(k)+dt;
u=A*cos(omega*t(k));
rho(:,:,k+1)=rho(:,:,k)+dt*(-1i)*((H+u*theta*mu)*rho(:,:,k)-rho(:,:,k)*(H+u*theta*mu));
rho11(k+1)=rho(1,1,k);
rho12(k+1)=rho(1,2,k);
rho21(k+1)=rho(2,1,k);
rho22(k+1)=rho(2,2,k);
y(k+1) = trace(P*rho(:,:,k+1));
% Solving for rho_hat %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
drho_hat_dt = -1i*((H+u*theta*mu)*rho_hat(:,:,k) - rho_hat(:,:,k)*(H+u*theta*mu))*GAMMA*(y(k)-y_hat(k))*(P*rho_hat(:,:,k) ...
+rho_hat(:,:,k)*P - 2*trace(P*rho_hat(:,:,k))*rho_hat(:,:,k));
rho_hat(:,:,k+1)=rho_hat(:,:,k)+dt*drho_hat_dt ;
% extracting matrix elements
rho_hat11(k+1)=rho_hat(1,1,k);
rho_hat12(k+1)=rho_hat(1,2,k);
rho_hat21(k+1)=rho_hat(2,1,k);
rho_hat22(k+1)=rho_hat(2,2,k);
% Solve for y_hat %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_hat(k+1) = trace(P*rho_hat(:,:,k));
% solving for theta_hat
dtheta_hat_dt = -1i*gamma*u*trace(P*(mu*rho_hat(:,:,k)-rho_hat(:,:,k)*mu))*(y(k)-y_hat(k));
theta_hat(k+1) = theta_hat(k) + dt*dtheta_hat_dt;
% Evaluating the error
e(k+1) = y(k) - y_hat(k);
end
subplot(223)
plot(t,abs(theta_hat));
grid on;
xlabel('Time [sec]','fontsize',20),
ylabel('$\hat{\theta}$','interpreter','latex','fontsize',20);
title("Plots",'fontsize',16);
subplot(224)
plot(t,abs(e));
grid on;
xlabel('Time [sec]','fontsize',20),
ylabel('$e$','interpreter','latex','fontsize',20);
title("Plots",'fontsize',16);
##### 0 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

### Respuestas (1)

Vaibhav el 12 de Abr. de 2024
Hi Kaydian
The issue with the convergence of "theta_hat" to 1 and e to 0 in your code may be due to several factors, including the learning rate "gamma", the initial conditions, the model used for the system and observer, or numerical stability.
Here are a few suggestions that might help improve convergence:
1. Check the Adaptation Gain (gamma): The value of "gamma" significantly affects the speed of convergence and stability. If "gamma" is too small, convergence will be slow, and if it's too large, the system might become unstable or oscillate.
2. Initial Conditions: The choice of initial conditions, especially for "theta_hat", can affect the convergence. You've already set theta_hat(1) = 1.5;, which is close to the true value of theta = 1. However, depending on the system dynamics and adaptation law, starting closer or exactly at the true value might not always be necessary or beneficial. Experiment with different initial conditions for "theta_hat" and "psi_hat".
3. Numerical Stability: The choice of "dt" (time step) can impact numerical stability. While your dt = 0.0001 seems small, ensuring numerical stability involves more than just the time step size; it also includes how the differential equations are solved. You might want to experiment with slightly larger "dt" values to see if it affects convergence without sacrificing too much precision.
Hope this helps!
##### 0 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

### Categorías

Más información sobre General Physics en Help Center y File Exchange.

### Community Treasure Hunt

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

Start Hunting!

Translated by