How to construct a function with given local minima? I need to construct a function $f(x,y)$ in which there are 3 minima: 2 local and 1 global as given below.

Locals are: z = f(0.2,0.3) = 0.7 | z = f(0.6,0.8) = 0.8
Global is: z = f(0.85,0.5) = 0.6

As a procedure, I've tried MLP, RBFNN and simple polynomial regression. But they just go through given x,y,z points without considering them as minima.
If you have any suggestions, it's more than welcome. Thanks in advance
 A: This thread may be useful to statisticians and data scientists because it shows how to construct arbitrarily "nasty" functions (that nevertheless are easy to handle mathematically and computationally) for testing algorithms that rely on optimization.

One way to construct functions with specified local properties (like local minima) is to assemble them from individual non-overlapping pieces.  A "piece" is a function that varies only within a specified region.
We may begin this construction by considering local functions of a single variable.  A nice candidate, due to its simplicity, ease of calculation, and control over its degree of differentiation, is proportional to a Beta density.  I will parameterize this to be supported on the interval $[-1,1],$ rising from zero to a maximum of $1$ at $0:$

$$\epsilon(x, a) = \left\{\eqalign{(1-x^2)^a & \quad |x| \lt 1 \\ 0 & \quad\text{otherwise.}}\right.$$

For $a\ge 1,$ this function has $2\lfloor a\rfloor-1$ continuous derivatives everywhere.

Its support can be scaled to the interval $[p-r,p+r]$ around $p$ simply by applying $\epsilon$ to $(x-p)/r.$  For $p, x\in\mathbb{R}^n,$ apply $\epsilon$ to the distance between $p$ and $x,$ written $|x-p|.$  Thus, the function

$$\psi_{p, r, a}(x) = \epsilon\left(\frac{|x-p|}{r}, a\right)$$

is nonzero only for points $x$ within distance $r$ of $p,$ attains a unique maximum of $1$ at $x=p,$ and its degree of differentiability is controlled by $a.$
Given a specification like that in the question--namely, a set of $d$ points $p_i$ and intended values $z_i$ at those points--a suitable linear combination of these functions will do the trick.  Merely find an $r \gt 0$ for which $2r$ is less than half the closest distance among those points, choose a suitable value of $a,$ and construct the function
$$f(x) = \sum_{i=1}^d z_i\, \psi_{p_i, r, a}(x).$$
If, in addition, you wish to guarantee that these points are local minima, take a linear combination of the negatives of the $\psi.$  Doing this requires specifying a global maximum $C$ strictly greater than any of the minima $z_i.$  Use

$$f(x) = C + \sum_{i=1}^d (z_i-C)\, \psi_{p_i, r, a}(x).$$

Why does this work?  First, at any point greater than distance $r$ from any $p_i,$ the sum vanishes and so $f$ is asymptotically constant at value $C.$  Second, the coefficients $z_i-C$ are all negative, demonstrating that (1) $C$ is the global maximum, (2) each $p_i$ is a local minimum, and (3) attains the value $C + (z_i-C) = z_i$ at $p_i$ (because the value of $\psi$ at $p_i$ is $1$).  Finally, there are no other local minima because each $\psi$ has a unique local maximum.
A similar method can be used to insert local maxima at specified points.

In the example in the question, here is a pseudo-3D plot of the result (with $a=2$ and $C=1$), looking in the positive $y$ direction from underneath:

Here is an image of the same function showing three cross-sections passing through the local minima:

This plot of the values of $f$ along the cross-sections (keyed by color) includes the three individual local minima shown as bubbles to demonstrate this truly is a solution to the problem:

In case any remaining details remain unclear, here is the R code that created these plots.  Although written for the two-dimensional case, it is readily modified for any number of dimensions (apart from the plots, of course).
#
# Define the basic "kernel" function.
#
eps_ <- function(x, a=2) (1 - pmin(1, abs(x))^2)^a
eps <- function(x, scale=1, a=2) eps_(x/scale, a)
psi <- function(x, scale=1, a=2) eps(sqrt(sum(x^2)), scale, a)
#
# Construct a linear combination of the kernels to achieve a specified set 
# of local minima as given by the matrix `X`.
#
X <- cbind(x = c(2/10, 6/10, 85/100),
           y = c(3/10, 8/10, 5/10),
           z = c(7/10, 8/10, 6/10))

f <- Vectorize(function(x, y, r=1, C=0) {
  C + sum(apply(X, 1, function(u) {
    psi(c(x,y) - u[1:2], scale=r) * (u[3] - C)
  }))
}, c("x", "y"))
#
# Check that the intended values are achieved at the points.
#
r <- min(dist(X[, 1:2]))/2
C <- max(X[, "z"]) + 0.2     # Choose a global maximum
print(mapply(f, X[, "x"], X[, "y"], r = r, C = C)) # Check
#
# Plot an image of the function.
#
x <- y <- seq(0,1,length.out=201)
z <- outer(x, y, f, r = r, C = C)
image(x, y, z, col=topo.colors(24))

