# Matlab code of Neural delay differential equation NDDE

10 views (last 30 days)
Commented: Muhammad on 5 Jan 2024
I have written a code of NDDE, but it is not correct and I am not able to simulate this code, because I am getting errors from it
my code
%% Dynamical System Modeling Using Neural ODE
% Parameters
par = [4; 2; 9.65];
tau = 1;
tau_max = 1.5;
num_trajectories = 5;
T = 30;
% Generating 100 trajectories
colors = jet(num_trajectories);
figure; % Create a single figure
for i = 0:(num_trajectories - 1)
c = 0.5 + i/99;
phi = @(t) c * ones(size(t));
g = @(t, y, Z) par(1) * Z / (1 + Z^par(3)) - par(2) * y;
% Solve the DDE using dde23 for baseline data
sol = dde23(@(t, y, Z) g(t, y, Z), tau, phi, [0 T]);
% Plot the current trajectory on the same figure
plot(sol.x, sol.y, 'Color', colors(i+1, :));
hold on;
end
title('Mackey-Glass Equation Trajectories');
xlabel('Time');
ylabel('x(t)');
hold off
numTimeSteps = 1000;
t = linspace(0, T, numTimeSteps);
X=deval(sol,t);
X_delay=interp1(t,X,t-tau);
xTrain = X(1,:)'; % Transpose to match your shape
% Plot x vs. delayed term
figure;
plot(X(1, :), X_delay);
xlabel('x(t)');
ylabel('x(t-\tau)');
title('x vs. delayed term');
grid on;
neuralOdeTimesteps = 10;
dt = t(2);
timesteps = (0:neuralOdeTimesteps)*dt;
nd=1;
tau = dlarray(tau_max*rand(1,nd));
nx = 1;
hiddenSize = 5;
% output_size = 1;
NDDE = struct;
NDDE.fc1 = struct;
sz = [hiddenSize nx*(1+nd)];
NDDE.fc1.Weights = initializeGlorot(sz);
NDDE.fc1.Bias = initializeZeros([sz(1) 1]);
NDDE.fc2 = struct;
sz = [hiddenSize hiddenSize];
NDDE.fc2.Weights = initializeGlorot(sz);
NDDE.fc2.Bias = initializeZeros([sz(1) 1]);
NDDE.fc3 = struct;
sz = [nx hiddenSize];
NDDE.fc3.Weights = initializeGlorot(sz);
tau = min(max(1e-5,tau),tau_max-1e-5);
learnRate = 0.002;
numIter = 1200;
miniBatchSize = 200;
plotFrequency = 10;
monitor = trainingProgressMonitor(Metrics="Loss",Info=["Iteration","LearnRate"],XLabel="Iteration");
Error using matlab.internal.lang.capability.Capability.require
This functionality is not available on remote platforms.

Error in matlab.ui.internal.uifigureImpl (line 32)
Capability.require(Capability.WebWindow);

Error in uifigure (line 34)
window = matlab.ui.internal.uifigureImpl(varargin{:});

Error in deepmonitor.internal.DLTMonitorView/createGUIComponents (line 159)
this.Figure = uifigure("Tag", "DEEPMONITOR_UIFIGURE");

Error in deepmonitor.internal.DLTMonitorView (line 116)
this.createGUIComponents();

Error in deepmonitor.internal.DLTMonitorFactory/createStandaloneView (line 8)
view = deepmonitor.internal.DLTMonitorView(model, this);

Error in deep.TrainingProgressMonitor/set.Visible (line 175)
this.View = this.Factory.createStandaloneView(this.Model);

Error in trainingProgressMonitor (line 140)
monitor.Visible = args.Visible;
numTrainingTimesteps = numTimeSteps;
trainingTimesteps = 1:numTrainingTimesteps;
plottingTimesteps = 2:numTimeSteps;
iteration = 0;
while iteration < numIter && ~monitor.Stop
iteration = iteration + 1;
% Create batch
[X, targets] = createMiniBatch(numTrainingTimesteps, neuralOdeTimesteps, miniBatchSize, xTrain);
% Evaluate network and compute loss and gradients
% Update network
% Plot loss
recordMetrics(monitor,iteration,Loss=loss);
% Plot predicted vs. real dynamics
if mod(iteration,plotFrequency) == 0 || iteration == 1
% Use ode45 to compute the solution
y=dde23(fun,par,delay,hist,time);
plot(xTrain(1,plottingTimesteps),xTrain(2,plottingTimesteps),"r--")
hold on
plot(y(1,:),y(2,:),"b-")
hold off
xlabel("x(1)")
ylabel("x(2)")
title("Predicted vs. Real Dynamics")
legend("Training Ground Truth", "Predicted")
drawnow
end
updateInfo(monitor,Iteration=iteration,LearnRate=learnRate);
monitor.Progress = 100*iteration/numIter;
end
function X = model(fun,par,delay,hist,time)
X = dde23(fun,par,delay,hist,time);
end
function dx=fun(t,y,Z,par)
dx = par.fc3.Weights*tanh(par.fc2.Weights*tanh(par.fc1.Weights*[y;Z]+par.fc1.Bias)+par.fc2.Bias);
end
% % Compute predictions.
X = model(fun,par,delay,hist,time);
% Compute L1 loss.
loss = l1loss(X,targets);
end
function [x0, targets] = createMiniBatch(numTimesteps,numTimesPerObs,miniBatchSize,X)
% Create batches of trajectories.
s = randperm(numTimesteps - numTimesPerObs, miniBatchSize);
x0 = dlarray(X(s,:));
targets = zeros([size(X,1) miniBatchSize numTimesPerObs]);
for i = 1:miniBatchSize
targets( 1:numTimesPerObs,i) = X( s(i) + 1:(s(i) + numTimesPerObs),:);
end
end
function bias = initializeZeros(sz)
bias = zeros(sz,'single');
bias = dlarray(bias);
end
function weights = initializeGlorot(sz)
Z = 2*rand(sz,'single') - 1;
bound = sqrt(6 / (sz(2)+ sz(1)));
weights = bound * Z;
weights = dlarray(weights);
end

Ben on 5 Jan 2024
I notice that the model function uses dde23. Unfortunately dde23 is not supported by dlarray and so you can't use this with automatic differentiation. The dlode45 solver is the only ODE solver we have that supports dlarray and automatic differentiation.
It may be possible to work something out with dlode45 but it could be quite complex. For example, for a delay differential equation with 1 delay and for , it's possible to reduce this to an ordinary differential equation for by writing for . Here's a proof of concept using .
p = 0.1;
tau = 1;
phi = @(t) 1;
F = @(t,y,z,p) -y + p*z;
% Solve with dde23 for fixed p on [0,tau]
sol = dde23(@(t,y,z) F(t,y,z,p), tau, phi, [0,tau]);
% Construct equivalent ODE on [0,tau]
G = @(t,y,p) F(t,y,phi(t),p);
% Solve for fixed p
odeSol = ode45(@(t,y) G(t,y,p), [0,tau], phi(0));
% Here I compute gradients of y(tau) with respect to params as a proof of
% concept.
y = dlode45(G,tspan,y0,params,DataFormat="CB");
end
To extend beyond would get a little bit complex - potentially you could interpolate the solution on to define for using interp1 (since it is supported by dlarray and automatic differentiation) and use this to solve for for using dlode45. However I expect this could get quite difficult to implement.
Thank you for nice explaination and I got a point where I am doing mistake