ML/NN as Function Evaluator for further Optimization (maximization) - Practical Example

Question

I am working on a production optimization problem; a very similar idea to what is described by Vegard Flovik How to use machine learning for production optimization. The following image, taken from the referred post, summarizes it very well:

First step is obvious, and I do have a model in the form of machine learning or neural networks model. How would I go about the second step? How can I use the trained model as the function evaluator for further multi-dimensional nonlinear optimization (e.g. maximization) either via Scipy, Bayesian Optimization etc.?

I cannot seem to find a practical example. Having a closed-form analytical function as the objective of an optimization problem is well-established. The article Optimization with SciPy and application ideas to machine learning by Tirthajyoti Sarkar gives a few examples using Scipy, & introduces packages that do optimizations with bound constrains and more. Yet examples are quite simple (a closed-form mathematical function) and he only glosses over the extension of such idea to use NN as the objective function, I'm quoting:

You are free to choose an analytical function, a deep learning network (perhaps as a regression model), or even a complicated simulation model, and throw them all together into the pit of optimization.

Any leads/hints/links are appreciated!

---

[Appendix]

In order to have a concrete example, let's imagine we have a dummy data set with a set of feature and a imaginary ProductionYield that is a nonlinear combination of input variables:

import numpy as np
import pandas as pd

df = pd.DataFrame(columns=['Pressure','Temprerature','Speed','ProductionYield'])

df['Pressure'] = np.random.randint(low= 2, high=10, size=2000)
df['Temprerature'] = np.random.randint(10, 30, size=2000)
df['Speed'] = np.random.weibull(2, size=2000)
df['ProductionYield'] = (df['Pressure'])**2 + df['Temprerature'] * df['Speed'] + 10
df['ProductionYield']= df['ProductionYield'].clip(0, 100)

   Pressure  Temprerature     Speed  ProductionYield
0         7            20  1.810557        95.211139
1         2            29  0.674221        33.552409
2         8            17  0.537533        83.138065
3         3            24  1.945914        65.701938
4         6            23  0.514679        57.837610

1.Predictive Algorithm (a simple Neural Network):

## Train/Test Split
from sklearn.model_selection import train_test_split

x_train, x_test, y_train, y_test = train_test_split(df[['Pressure','Temprerature','Speed']].values, df['ProductionYield'].values, test_size=0.33, random_state=42)

## Build NN Model
import tensorflow as tf
from tensorflow.keras import layers

def build_model():
    
    # create model
    model = tf.keras.Sequential()
    model.add(layers.Dense(64, input_dim=3, kernel_initializer='normal', activation='relu'))
    model.add(layers.Dense(128, kernel_initializer='normal', activation='relu'))
    model.add(layers.Dense(1, kernel_initializer='normal'))
    
    # Compile model
    model.compile(loss='mean_squared_error', optimizer='adam')
    
    return model

model = build_model()
model.fit(x_train, y_train,
          validation_split=0.2,
          verbose=0, epochs=1000)

2.Otimization [Core of the Problem]:

Problem lies herein, when a ML/NN is trained, I do not get to see (export as I wish) the mathematical form of the function (here in this example NN) and its variables (which should be my feature variables) to do the optimization as we do with closed-form explicit mathematical functions.

[UPDATE 15.01.2021]:

Following Valentin's great answer, I've put pieces together in a practical example showcasing how one can use a ML/NN model as an input function for further optimization (herein via scipy.optimize) using the dummy data set shown in the Appendix. Please see this notebook for more details.

Answer 1

This post seems similar to yours and may help. It seems that what you are looking for is a derivative-free optimization method. The Wikipedia page for the concept lists such methods.

Intuitively, these techniques will sample the function (the network in your case) with various inputs (pressure, temperature, speed) and will figure out which inputs optimize it. Where they differ is in their sampling strategy, as it may be impractical or expensive to sample.

You can use scipy.optimize.minimize to do that. Pass your network as func and use an initial guess, which can be the last values of your variables. Scipy expects a function with the following signature: fun(x, *args) -> float where x is a one-dimensional numpy array. This might mean you will need to wrap your network in something like this:

def wrapper(x, *args) -> float:
    
    network_input = _numpy_to_valid_network_input(x)

    network_output = network.predict(network_input, *args)

    scipy_output = _network_output_to_float(network_output)

    return scipy_output

And then, you can pass wrapper as your func. Negating the output of scipy_output will turn the minimization problem into a maximization problem. If your input variables are bounded, i.e. [0, 100], you can do two things:

• Use algorithms that allow you to define these bounds explicitly (i.e. L-BFGS-B using the bounds parameter)
• Implicitly bound your inputs by using a sigmoid function for instance. By doing something like I did below, you can create a function bounded_value that always returns values within the range of your choice, even if the optimization algorithm might try any float.

import numpy as np

def sigmoid(x):
  return 1 / (1 + np.exp(-x))

def bounded_value(x, min_value=0, max_value=100):
   return min_value + (sigmoid(x) * (max_value - min_value))

If you want to bound your output, that's even easier. If your goal is to minimize, return a very large value if your network outputs any value outside of your range. Obviously, ensuring that your network always returns sensible values could be achieved by tweaking the loss used during training (which is difficult with AutoKeras).

If you want to implement your own method, I would recommend coordinate descent as you don't have many input dimensions (3) and it is quite simple to implement from scratch. Obviously, a brute-force approach where you randomly sample the space and choose the inputs that yielded the best function value is even easier to implement.