Constructing function - f(x,y) for the given minimums (Python) Problem Statement:
I need to construct a function f(x,y) in which there're 3 minimums. 2 local and 1 global which are written 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


What I've tried so far
I've tried Polynomial regression to construct features that estimate the given minimum points. However, an estimated function z = f(x,y) just going through all the given minimum points without considering them as minimums. I need an approach by means of which I can construct a function that has 3 minimums that are mentioned above. I'd be more than happy if you could provide a code snippet of your approach. As a programming language, I use Python.
Code
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

degree = 6


def func(x1, x2):
    arr = []
    for i in range(degree, 0, -1):
        arr.append(x1 ** i)

    for i in range(degree, 0, -1):
        arr.append(x2 ** i)

    arr.append(1)

    return np.array(arr)


def der_func(x1, x2):
    arr1 = []
    arr2 = []
    for i in range(degree - 1, -1, -1):  
        arr1.append((i + 1) * x1 ** i)  

    for i in range(degree + 1):
        arr1.append(0)

    for i in range(degree):
        arr2.append(0)

    for i in range(degree - 1, -1, -1):
        arr2.append((i + 1) * x2 ** i)

    arr2.append(0)

    return np.array([arr1, arr2])


given_x = np.array([[0.2, 0.3],
                    [0.6, 0.8],
                    [0.85, 0.5]])

given_y = np.array([0.7, 0.8, 0.6])

data_x = []
data_y = []

for i in range(len(given_x)):
    x1x2 = given_x[i]

    data_x.append(func(x1x2[0], x1x2[1]))
    data_y.append(given_y[i])

    if i &lt; 3: # in case there's more than 3 data point in given_x
        data_x.append(der_func(x1x2[0], x1x2[1])[0])
        data_y.append(0)

        data_x.append(der_func(x1x2[0], x1x2[1])[1])
        data_y.append(0)

data_x = np.array(data_x)
data_y = np.array(data_y)

w = np.linalg.inv(data_x.T @ data_x) @ data_x.T @ data_y
pred_y = data_x @ w

cost = np.mean((data_y - pred_y) ** 2)
print(&quot;Cost | &quot;, cost)

#######################################

X = np.linspace(0, 1, 100)
Y = np.linspace(0, 1, 100)

X, Y = np.meshgrid(X, Y)
Z = np.zeros(np.shape(X))
for i in range(len(X)):
    for j in range(len(X[0])):
        Z[i, j] = np.array([func(X[i, j], Y[i, j])]) @ w

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.contour(X, Y, Z, 150)
ax.scatter(given_x[:3, 0], given_x[:3, 1], given_y[:3])
plt.show()
One way to go is to define functions that are non-zero only near your 'minimal' points, and add them. The goal is to avoid overlap, such that a function associated with a point won't modify the value of your global function around other points. 

import numpy as np

def func0(x1):
    lim = 0.01
    dist = x1*x1
    if dist&lt;lim:
        value = np.exp(-dist/(lim-dist))
    else: 
        value = 0
    return value


That will give you a 1D function with nice properties (value 1 at 0, 0 outside the neigborhood of 0) : 

x = np.linspace(-1, 1, 1000)
z = np.array([func0(i) for i in x])
Z = z.reshape(1000)

plt.plot(Z)
plt.show()




Modify it to get a 2 D function centered in (a,b) :

def func(x1,x2):
    a = 0
    b = 0
    lim = 0.01
    dist = (x1-a)*(x1-a) + (x2-b)*(x2-b)
    if dist&lt;lim:
        value = np.exp(-dist/(lim-dist))
    else: 
        value = 0
    return value 


This will give you a 2d function :

x = y = np.linspace(-1, 1, 100)
z = np.array([func(i,j) for j in y for i in x])
Z = z.reshape(100, 100)

x, x = np.meshgrid(x, y)

plt.imshow(Z,extent=(x.min(), x.max(), y.min(), y.max()))
plt.colorbar()
plt.show()




Then build the threes functions associated with your points. 

def func1(x1,x2):
    a = 0.2
    b = 0.3
    lim = 0.01
    dist = (x1-a)*(x1-a) + (x2-b)*(x2-b)
    if dist&lt;lim:
        value = np.exp(-dist/(lim-dist))
    else: 
        value = 0
    return value

def func2(x1,x2):
    a = 0.6
    b = 0.8
    lim = 0.01
    dist = (x1-a)*(x1-a) + (x2-b)*(x2-b)
    if dist&lt;lim:
        value = np.exp(-dist/(lim-dist))
    else: 
        value = 0
    return value

def func3(x1,x2):
    a = 0.85
    b = 0.5
    lim = 0.01
    dist = (x1-a)*(x1-a) + (x2-b)*(x2-b)
    if dist&lt;lim:
        value = np.exp(-dist/(lim-dist))
    else: 
        value = 0
    return value


And add them (up to coefficients that determine the value on points of interest) to form a global function : 

def global_func(x,y):
    return 1-0.3*func1(x,y)-0.2*func2(x,y)-0.4*func3(x,y)


This will give you :

x = y = np.linspace(0, 1, 1000)
z = np.array([global_func(i,j) for j in y for i in x])
Z = z.reshape(1000, 1000)

x, x = np.meshgrid(x, y)

plt.imshow(Z,extent=(x.min(), x.max(), y.min(), y.max()))
plt.colorbar()
plt.show()




Finally you can check that :

print(global_func(0.2,0.3))
print(global_func(0.6,0.8))
print(global_func(0.85,0.5))


returns :

0.7
0.8
0.6


You have a continuous function with all expected properties by construction. 

If you are unhappy with it being constant on a lot of the domain, you can modify the individual functions so that they don't stop at the limit with :

value = np.exp(-dist)


And in your global function, determine which 'minimal' point is the nearest and only apply a function of the form 1-coef * funci based on which point is nearest. This way the individual functions won't overlap and thus won't modify value of the global function around other 'minimal' points.