Azzera filtri
Azzera filtri

How do I get this code to converge?

1 visualizzazione (ultimi 30 giorni)
Kaydian Quintero
Kaydian Quintero il 20 Lug 2021
Risposto: Vaibhav il 12 Apr 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);

Risposte (1)

Vaibhav
Vaibhav il 12 Apr 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!

Categorie

Scopri di più su General Physics in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by