How do I get this code to converge?

1 次查看（过去 30 天）
Kaydian Quintero
Kaydian Quintero 2021-7-20
回答： Vaibhav 2024-4-12
% 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);

请先登录，再回答此问题。

回答（1 个）

Vaibhav
Vaibhav 2024-4-12
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!

类别

在 Help Center 和 File Exchange 中查找有关 ThingSpeak 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by