Constructing function - f(x,y) for the given minimums (Python)
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.