This is a perfect problem for active learning. Methods based on Bayesian Optimization are particularly powerful for optimizing black-box functions which are expensive to evaluate (i.e. running an experiment in the lab). There are a few BO packages out there which may be of interest, Martin Kraisser's blog has a nice overview.
I noticed that the features in your last experiment don't add up to 1 which I am assuming was a typo. For the demo I changed that entry to x2 = 0.6.
Here is a sample I threw together in python using GPyOpt, a Gaussian Process based package:
import numpy as np
import GPyOpt
x_init = np.array([[0.9,0.0,0.1,0.0],
[0.0,0.9,0.1,0.0],
[0.45,0.45,0.0,0.1],
[0.6,0.3,0.05,0.05],
[0.3,0.6,0.05,0.05]])
y_init = np.array([[16.5],[8.6],[12.6],[18.9],[9.8]])*(-1)
domain = [{'name': 'x1', 'type': 'continuous', 'domain': (0,1.0)},
{'name': 'x2', 'type': 'continuous', 'domain': (0,1.0)},
{'name': 'x3', 'type': 'continuous', 'domain': (0,1.0)},
{'name': 'x4', 'type': 'continuous', 'domain': (0,1.0)}
]
constraints = [
{'name':'const_1', 'constraint': '(x[:,0] + x[:,1] + x[:,2] + x[:,3]) - 1 - 0.001'},
{'name':'const_2', 'constraint': '1 - (x[:,0] + x[:,1] + x[:,2] + x[:,3]) - 0.001'}
]
bo_step = GPyOpt.methods.BayesianOptimization(
f = None,
domain = domain,
constraints = constraints,
X = x_init,
Y = y_init
)
x_next = bo_step.suggest_next_locations()
print(x_next)
print(np.sum(x_next))
Notes:
(1) GPyOpt only accepts constraints in a certain form, that's why there are 2 which constrain y on the interval [0.999,1.001].
(2) GPyOpt is set up to minimize functions. Therefore I multiplied your target values by -1.
This example suggests that your next experiment should be run at:
x1 = 0.69
x2 = 0.21
x3 = 0.08
x4 = 0.02
BO algorithms can be tuned to give different results based on your preferences for exploiting existing information versus exploring new areas of the space. I'm not sure what GPyOpts standard settings are so If you are interested it could be worth looking at the documentation.