Matlab code of Neural delay differential equation NDDE

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);
gradDecay = 0.9;
sqGradDecay = 0.999;
learnRate = 0.002;
numIter = 1200;
miniBatchSize = 200;
plotFrequency = 10;
averageGrad_p = [];
averageSqGrad_p = [];
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
    [loss,grad_par,grad_tau] = dlfeval(@modelLoss,timesteps,X,par,targets);
    % Update network
    [par,averageGrad_p,averageSqGrad_p] = adamupdate(par,grad_par,averageGrad_p,averageSqGrad_p,iteration, learnRate,gradDecay,sqGradDecay);
    [tau,averageGrad_t,averageSqGrad_t] = adamupdate(tau,grad_tau,averageGrad_t,averageSqGrad_t,iteration, learnRate,gradDecay,sqGradDecay);
    % 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
function [loss,grad_par,grad_tau] = modelLoss(tspan,X0,par,targets)
% % Compute predictions.
X = model(fun,par,delay,hist,time);
% Compute L1 loss.
loss = l1loss(X,targets);
% Compute gradients.
grad_par = dlgradient(loss,par);
grad_tau = dlgradient(loss,tau);
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

Answer 1

[Ben]

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));

% Compute gradients with dlode45

% Here I compute gradients of y(tau) with respect to params as a proof of

% concept.

[y,grad] = dlfeval(@solveAndGradient,G,[0,tau],dlarray(phi(0)),dlarray(p));

function [y,grad] = solveAndGradient(G,tspan,y0,params)

y = dlode45(G,tspan,y0,params,DataFormat="CB");

grad = dlgradient(y,params);

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.

[Muhammad]

Thank you for nice explaination and I got a point where I am doing mistake