lululxvi / deepxde

A library for scientific machine learning and physics-informed learning
https://deepxde.readthedocs.io
GNU Lesser General Public License v2.1
2.64k stars 739 forks source link

Use `auxilary_vars` at prediction to make new predictions #443

Closed NouamaneTazi closed 2 years ago

NouamaneTazi commented 2 years ago

Hello @lululxvi , thank you for the amazing library! I would like to solve an optimal control problem.

v(t)=R.i(t) + dΦ(t)/dt

where v(t) is the controller.

I trained an initial model inspired from Lorenz inverse forced example, which fitted very well to the data

R = dde.Variable(1.0)

interpolate time / lift vectors (for using exogenous variable without fixed time stamps)

def ex_func2(t): spline = sp.interpolate.Rbf( ex_time, ex_input, function="thin_plate", smooth=0, episilon=0 ) return spline(t[:, 0:])

def system(t, y, ex): """system. v = ex (exogenous) v = R.i + dphi/dt """ i, phi = y[:, 0:1], y[:, 1:2] dphi_t = dde.grad.jacobian(y, t, i=1) return ex - R*i - dphi_t

def boundary(_, on_initial): return on_initial

geom = dde.geometry.TimeDomain(0, 0.2)

Initial conditions

ic1 = dde.IC(geom, lambda X: ob_y[0,0], boundary, component=0) ic2 = dde.IC(geom, lambda X: ob_y[0,1], boundary, component=1)

Get the train data

observe_i = dde.PointSetBC(observe_t, ob_y[:, 0:1], component=0) observe_phi = dde.PointSetBC(observe_t, ob_y[:, 1:2], component=1)

data = dde.data.PDE( geom, system, [ic1, ic2, observe_i, observe_phi], num_domain=1000, num_boundary=2, anchors=observe_t, auxiliary_var_function=ex_func2, )

net = dde.maps.FNN([1] + [40] * 3 + [2], "tanh", "Glorot uniform") model = dde.Model(data, net) model.compile("adam", lr=0.001, external_trainable_variables=[R])

variable = dde.callbacks.VariableValue(R, period=1000) losshistory, train_state = model.train(epochs=60000, callbacks=[variable])

* Inference:
```python
t_true, y_true = get_data(filename)
y_pred = model.predict(t_true)

# Plot results
plt.plot(t_true, y_true[:,1], label="flux")
plt.plot(t_true, y_pred[:,1], label="pred")
plt.legend()
plt.figure()
plt.show()

image

I would like now to use different values for v(t) and have new predictions. Any ideas how to do that? Thanks.

NouamaneTazi commented 2 years ago

I tried modifying model.data.auxiliary_var_fn to use a new function and retrain the model for a couple epochs. But it doesn't help the training at all. :/

filename='phase2/input0.csv'

ex_time, ex_input = get_ex_data(filename) # exogenous data

# interpolate time / lift vectors (for using exogenous variable without fixed time stamps)
def ex_func3(t):
    spline = sp.interpolate.Rbf(
        ex_time, ex_input, function="thin_plate", smooth=0, episilon=0
    )
    return spline(t[:, 0:])
data.auxiliary_var_fn = ex_func3 # I did check that model.data.auxiliary_var_fn was updated correctly

losshistory, train_state = model.train(epochs=10000, callbacks=[variable])

image

lululxvi commented 2 years ago

You need to re-do everything for new v(t), unless you would like to do something like transfer learning

NouamaneTazi commented 2 years ago

For now I managed to "rewire" the voltage so that I give it to the NN as an input. But I don't think the geometry is intended to work this way (I basically only train on the anchor points, completely ignoring the space geometry). But I'm very doubtful about the loss I get in this case (although the PDE residuals are far from 0, I get a training loss (the first one near 1e-6)

R = dde.Variable(1.0)

def system(inputs, y):
    """system.
    v is second input
    v = R.i + dphi/dt
    """
    t, v = inputs[:, 0:1],inputs[:, 1:2]
    i, phi = y[:, 0:1], y[:, 1:2]
    dphi_t = dde.grad.jacobian(y, inputs, i=1, j=0)
    return v - R*i - dphi_t

def boundary(_, on_initial):
    return on_initial

geom = dde.geometry.Disk(0, 1)
timedomain = dde.geometry.TimeDomain(0, 0.2)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

# Initial conditions
ic1 = dde.IC(geom, lambda X: ob_y[0,0], boundary, component=0)
ic2 = dde.IC(geom, lambda X: ob_y[0,1], boundary, component=1)

# Get the train data
observe_i = dde.PointSetBC(observe_inp, ob_y[:, 0:1], component=0)
observe_phi = dde.PointSetBC(observe_inp, ob_y[:, 1:2], component=1)

data = dde.data.TimePDE(
    geomtime,
    system,
    [observe_i, observe_phi],
    num_domain=0,
    num_test=100,
    anchors=observe_inp,
)

# net = dde.maps.FNN([3] + [40] * 5 + [2], "tanh", "Glorot uniform")
net = dde.maps.ResNet(3,2,50,1, "tanh", "Glorot uniform")
model = dde.Model(data, net)
model.compile("adam", lr=0.001, external_trainable_variables=[R])

image

lululxvi commented 2 years ago

It is OK to reuse geometry in that way. Is your result good?