Using the deepxde.data.constraint function

[9erik1]

Hello, I was just wondering if anyone can provide a basic example of how one might use the data.constraint function of deepxde (see https://deepxde.readthedocs.io/en/latest/_modules/deepxde/data/constraint.html?highlight=constraint). From the documentation, all I see is a link to the source code, which I don't understand because I'm not well-versed enough in python and defining classes. I'm presuming this is what's needed to place constraints on the data generated by the program? I've looked for similar questions but have only found specific answers that don't seem to pertain to what I'm doing.

If it helps to give the background, I'm working with a group on trying to solve the equation of motion for the two body problem, written as the effective one body problem. The equations of motion, in dimensionless form, is:

• (dr/dt)^2 = [(e^2 - 1) + 1/r - 1/r^2] = (e*sin(theta))^2 -- and --
• dtheta/dt = 1/2r^2

where (r, theta, e) are the radius, polar angle and eccentricity. If I use (dr/dt)^2 = [(e^2 - 1) + 1/r - 1/r^2] as my loss, then I get bad results -- mainly, the algorithm always runs to a circular orbit. If I use dr/dt = e*sin(theta), then the algorithm seems to (a) prefer negative radii/theta and (b) not want theta to change signs, and prefers a solution where theta increases asymptotically towards zero. My guess is that the squared terms are causing problems, as well as the fact that the diff eqn for radius yields a minimum radius, which the algorithm is trying desperately to avoid going below.

I put a picture below to illustrate -- I have the initial theta set to -pi, and as you can see theta increases asymptotically towards zero while r seems to exponentially decay to compensate. Additionally, if I set the initial condition on theta to be somewhere were sin(theta) >= 0, then the solution tends to jump immediately to a negative theta and then again increase asymptotically towards zero.

Thanks for your help!

[lululxvi]

Your problem is an ODE system, so you can use dde.data.PDE.

[9erik1]

Hi Prof. Lu, thanks for the response -- so we are currently using dde.data.PDE (I've pasted my code below if it helps), essentially just testing an altered copy of the example given in https://deepxde.readthedocs.io/en/latest/demos/pinn_forward/ode.system.html. What I think is the problem, is that the ODE has solutions that are unphysical, but the network doesn't realize this. For the orbit equation, there is a point of closest approach (in the case of an elliptical orbit, this would be at theta = 0), and it seems to have a hard time providing a solution around this point. In some formulations of the problem, I get half an ellipse (see picture below), which reflects the radius/angle solution diagram I put in my original post.

Basically, we're dealing with a DE that potentially has (dr_dt)^2 < 0, so we want to avoid that by constraining the solutions to the problem. I can't figure out how to do that with dde.data.PDE (if it is even possible), and was assuming that dde.data.constraint may be the way to go.

def ode_system(t, S):
    """
    t = time (dimensionless); 
    S = state (dimensionless radius and angle)
    """
    
    r, th = S[:,0:1], S[:,1:] # dimensionless units, but by design 
    #tf.print(r)
    r_dot = dde.grad.jacobian(S, t, i=0) 
    th_dot = dde.grad.jacobian(S, t, i=1) # angular velocity
    r_dot_loss = r_dot - eps*tf.math.sin(th) # r_dot = eps*sin(th)
    th_dot_loss = 2*th_dot*r**2 - 1 # th_dot = 1/2/r**2
    return [r_dot_loss, th_dot_loss]

def boundary(_, on_initial):
    return on_initial

geom = dde.geometry.TimeDomain(0, 500)
ic2 = dde.icbc.IC(geom, lambda x: -1*pi, boundary, component=1)
data = dde.data.PDE(geom, ode_system, [ic2], 100, 3, solution=None, num_test=1000)

layer_size = [1] + [50] * 3 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.FNN(layer_size, activation, initializer)

model = dde.Model(data, net)
print(model)
model.compile("adam", lr=0.001)#, metrics=["l2 relative error"])
losshistory, train_state = model.train(iterations=10000)

dde.saveplot(losshistory, train_state, issave=True, isplot=True)

[lululxvi]

I don't fully understand your question.

[9erik1]

Hello Prof. Lu,

So this is an orbit equation, and setting the eccentricity e should determine the shape of the orbit. In my screenshots above, I think I used e = 0.9 (highly elliptical orbit). We can't recover the full proper solution with the neural network (the solutions change depending on how I set up the loss terms, but for the most straightforward ways the network seems to panic and settle on a circular orbit). I believe it has to do with the facts that:

(a) the equation of motion deals with distance r. We know that r > 0, but the network doesn't. Sometimes it spits out solutions with negative r(t) function.
(b) there is a point of closest approach in any given orbit (set dr_dt = 0 in the first equation of motion below and solve quadratic equation). The r < 0 solutions of the quadratic equation must be rejected, but the network has no knowledge of this. The network seems to do poorly around the quadratic equation roots and wants to avoid them.

The point is that the equations of motion only work with a-priori knowledge that constrains the solutions. The neural network has no such knowledge, and I think this is the problem. My thinking is that if we can somehow constrain the solutions beforehand, we may be able to recover the proper solution. I've tried adding a term to the normal losses that will be zero if r > r_closest > 0, but will be a really big number if r < r_closest, but this did not work.

In other words, I want to constrain the solution of r(t) so that it doesn't violate the above constraints, but I don't know the proper way to do this (or if there even is one) -- I tried searching the documentation and the best thing I could find was the dde.data.constraint page, but I don't understand how this is supposed to work.

Equations of motion:

• (dr_dt)**2 = (e**2 - 1) + 1/r - 1/4/r**2 (sorry there was a slight typo in my original post)
-- and --
• dtheta_dt = 1/2/r**2

[lululxvi]

If your problem is PDE, then you need to use dde.data.PDE or dde.data.TimePDE. Extra constraints can be implemented, but the implementation depends on your constraint. For r>0, you can check the example at https://epubs.siam.org/doi/10.1137/21M1397908 and the code there.

[cnn-yangdx]

Hello Prof. Lu, when runing deepxde for a inversion problem, how to constain values of external_trainable_variables, and let them to be positive during training, thank you so much