Optimal Predictor Desgin

This repository contains codes and results for optimal estimation of heat equation by means of shape optimization and neural networks.

Consider a one-dimensional heat-conductive bar over the interval . Let be the temperature of the bar at location and time . The changes in the temperature is governed by the equation:

Forward simulation

Forward simulation involves a function that gives the output to the model given the inputs. For our specific example, the inputs are an initial temperature profile and a sensor shape set . The output is the temperature measured by a sensor in the set ; that is, .

There are various methods to solve the heat equation and find the solution for every initial condition. We use forward-time central-space finite-difference discretization method to find the solution of the heat equation. The following Python function is created that yields the solution

u = FTCS(dt,dx,t_max,x_max,k,u0)

The parameters of the function are defined below.

Time increment:	Space discretization:	Final time:	Length of the bar:	conductivity:
0.1	0.01	100	1	0.0003

For the specified parameters and the following initial condition , the solution is obtained by running the code

>> .\forward_sim.py

Neural-Network Estimator

A neural-network estimator is trained from some set of initial conditions to estimate the solution of the heat equation for any arbitrary initial condition. The set of initial conditions selected for training is

where is changed from 1 to N to create new training sample. The training inputs are illustrated in the following figure.

Use the function generate_data() to create training data. Training data are also stored in the external files input.npy and output.npy.

Building a model

We build a model using keras library of tensorFlow. The estimator is indicated by model and is consitruced in steps as follows.

1. Choosing Layers

Sequential Layer

A sequential layer is defined using the comand tensorflow.keras.Sequential. Layers are added one by one using the command model.add. Three layers are often present: Input Layer, Dense Layer, Output Layer.

model = Sequential()
model.add(Dense(100, input_dim = c1, activation='selu'))
model.add(Dense(500, activation='selu'))
model.add(Dense(1000, activation='selu'))
model.add(Dense(500, activation='selu'))
model.add(Dense(c, activation='selu'))

The architecture of this model is as follows

RNN Layer

An RNN layer is a recurrent neural network in which the output is fed back to the network. A schematic of the network is depicted below

RNN schematic	current architecture

To add a recurrent layer to a sequential model in keras, the layer SimpleRNN is used.

model = Sequential([
    SimpleRNN(c, activation=act_function, return_sequences=True, input_shape=[None, cl])
    ])

The following is the simulation result for an RNN predictor

CNN Layer

A CNN layer applies various filters to a time series and yields a time series width shorter with depending on the size of its filter. A schematic of the network is depicted below

CNN schematic	current architecture

To add a 1D convolutional layer to a sequential model in keras, Conv1D is used as follows

model = Sequential([
    Conv1D(filters=c, activation=act_function, kernel_size=cl,
     strides=1, padding="same", input_shape=(1, cl))
     ])

The following is the simulation result for a CNN predictor

2. Choosing the Compiler

Optimization method, loss function, and performance metrics are chosen in this step.

model.compile(optimizer='adam', loss='mean_squared_error', metrics=['accuracy'])

3. Training the Model

Sample data are fed to the model to train it. The data are divided into epochs and each epoch is further divided into batches.

model.fit(input, output, epochs=1, batch_size=m)

4. Evaluating the Performance

We use some test data to evaluate the performance of the trained model. This includes feeding some sample input and output to the model and calculate the loss and performance metric.

eval_result = model.evaluate(u0, u_real.reshape((1,c,r)), batch_size=1)
print('evaluation result, [loss, accuracy]=' , eval_result)

5. Estimation (i.e. Making Predictions)

For the time being, we assume and . Estimation is

u_pred=model.predict(np.asarray(u0).reshape((1,c)), batch_size=1)

An activation function can be chosen in each layer. An activation function should reflect the original mathematical model. Optimizer and loss fuctions are to be chosen as well.

Different initial conditions are tested to observe the performance of the estimator. We have cosidered the follwing initial conditions:

Shape Optimization

Let indicate the sensor shape, and be the real solution to the heat equation and be the prediction. Consider the following linear-quadratic cost function:

where is a characteristic function indicating the sensor region.

The derivative of with respect to is

The function is the derivative of with respect parameter . This derivative is

The function is the gradient of the network with respect to its input, and is the gradient of the network with respect to sensor locations. Also, the function is the derivative of with respect to sensor location . As a result, the gradient of the cost function with respect to sensor location is

We are interested in the discrete implementation of sensor shape optimization. In this implementation, we interpret as a vector with as many entries as the vector . Each entry is the probability of sensor presence at each node. As an example, shows that the probability of sensor presence at the first node is 0.5, at the second node is 0.1, and so on. The characteristic function is then the diag operation on omega. So, the term is replaced with the matrix multiplication diag(omega)*u_real. Also, the derivative is the derivative with respect to each probability of sensor presence. For instance, the derivative with respect to probability of sensor at node 1 is diag([1,0,0,...,0]). It is assumed that the gradient of the newrok with respect to omega is negligible.

shape_optimizer.py

Use the python function LQ_cost(u_pred,model,omega) to calculate the cost.

def LQ_cost(omega):

    omega_ = [i for i, value in enumerate(omega) if value < 0.5]

    input = array(training_data[0])
    input[:,omega_] = 0
    output = array(training_data[1])

    model.fit(input, output, batch_size=1000, epochs=4,verbose=0)

    u_pred = model.predict(matmul(u_real,diag(omega)))
    cost = (u_pred - u_real)**2
    cost = trapz(trapz(cost, dx=dx), dx=dt) 
    cost += 350*sum(omega) #(c-len(omega_))
    
    return cost

Keras uses the method keras.backend.gradients for calculating the gradient. Use the python function LQ_grad(omega) to calculate the gradient

def LQ_grad(omega):

    u_pred = model.predict(matmul(u_real,diag(omega)))

    grads = K.gradients(model.output, model.input)[0]
    gradient = K.function(model.input, grads)
    g = gradient(matmul(u_real,diag(omega)))

    D = 2*g*(u_pred-u_real)
     
    grad = zeros(c)
    for i in range(0,c):
        int = D[:,i]
        grad[i] = trapz(int, dx=dt) + 350
    
    return grad

The optimal actuator shape is then obtained by running the code

res = minimize(LQ_cost,omega_0,method='trust-constr',
           jac=LQ_grad,
           bounds=bounds,
           options={'verbose': 1, 'maxiter': 100, 'disp': True},
           callback=callback_cost)

Use the function animate(u_real,model,omega,FileName) to illustrate the performance of the estimator and save the result.

Time evolution for optimal sensor arrangement	Cost vs iterations

Time evolution when	Time evolution when

Time evolution when is random	Time evolution when is zero array

References

[1] M.S. Edalatzadeh and R. Herzog, "Optimal Prediction using Learning and Shape Optimization", arXiv:2105.05404.