Finding the optimal combination of inputs which return maximal output for a black box model

One of the challenges that I have been facing when applying ANN to regression tasks on my job is that in order to find the optimal out for a given range of inputs, I have to feed a multidimensional meshgrid to my model and then simply pick the highest value. However, this is overall a very computationally costly solution. The lenght of the text bellow might be scary but it just my attempt to better explain it.

Let me explain with other words. Supposing that I have 9 inputs for my ANN, and then I want to check which combinations of values of my features that returns me the highest outcome. I am currently overcoming the problem by just creating a 9D-mesh and simply predict the value for each sample and then identifying the optimal row. Nevertheless, this takes an exhaustive amount of time to work. Therefore, I am looking for a way be able to more efficiently reach this optimal output value, if possible at all.

In code, it would look something like this: (just a simple and made up example not really realistic in python):

import numpy as np
from itertools import product
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split 
import tensorflow.keras
from keras.models import Sequential
from keras.layers import Dense
import pandas as pd
import math

x1 = np.linspace(0,20,6)
x2 = np.linspace(0,20,6)
X = pd.DataFrame((product(*[x1,x2])))
y1 = 5*np.cos(np.deg2rad(X[0]))
y2 = 5 - 1*np.exp((-X[0]**2/np.deg2rad(10)**2)*np.cos(np.deg2rad(X[1])))
y = np.array([y1 + y2]).T

Setting a blackbox model, in this case, a neural network

x_scaler = MinMaxScaler()
y_scaler = MinMaxScaler()
X_scaled = x_scaler.fit_transform(X)
y_scaled = y_scaler.fit_transform(y)
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_scaled, test_size=0.2, random_state=0)

model = Sequential()
model.add(Dense(100, input_dim = 2, activation = 'relu'))
model.add(Dense(100, activation='relu'))
model.add(Dense(100, activation='relu'))
model.add(Dense(1, activation = 'relu'))
model.compile(optimizer = 'Adam', loss = 'mean_squared_error')
epochs_hist = model.fit(X_train, y_train, epochs = 100, batch_size = 50, validation_split = 0.2)

Now that I fit my model, I will use the meshgrid for several intervals in order to find the optimal in the specified range:

x1_pred = np.linspace(0,20,21) 
x2_pred = np.linspace(0,20,21)
X_pred = pd.DataFrame((product(*[x1_pred,x2_pred])))
X_pred_test = x_scaler.fit_transform(X_pred)
y_pred = model.predict(X_pred_test)
y_pred = y_scaler.inverse_transform(y_pred)

So, supposing that I doing something similar for reaching the optimal, but in this case with 9 inputs, then it is clear how computionally unfeasiable that calculation will be. Hence, it comes my question of how to find the optimal combination of inputs which return maximal output of a blackbox model such as ANN.

from Finding the optimal combination of inputs which return maximal output for a black box model