Advent Challenge 7: Balance the storage carts!
Python 3, 533 494 476 bytes

from numpy import*
def o(a,v,t=[],n=1):
 if len(v)&lt;1:return[[sum((d-e)**2for d,e in zip([sum(a)/sum(r)for a in zip(*[[a*k for k in b]for a,b in zip(r,t)])],[-~a/2 for a in a.shape]))**.5,a]]
 s,k=v[0],[[]]
 for p,i,j,l in[[copy(a),i,j,(i-~s/2,j-~s/2)]for i in range(a.shape[0]-s+1)for j in range(a.shape[1]-s+1)if all(a[i:i+s,j:j+s]&lt;1)]:p[i:i+s,j:j+s]=n;k+=o(p,v[1:],t+[l],n+1)
 return k
w=input()
q,r=w[0],w[1:]
print(min([a for a in o(zeros(q),r)if a],key=lambda a:a[0])[1])

Try it online!
I used the numpy library in Python.  It is a very popular Python library and it specializes matrix/multi-dimensional array manipulation.
Mine finds a different, but equally valid, solution.  I did double check my work and, when I tried to find all solutions, it did find the one presented in the question as well.
Basic methodology is to bruteforce it.  So, every possible present location is checked, then the one with the least distance between it's center of gravity and the center of the slay is returned (or at least one of the possible solutions).

Pseudo-Ungolfed Version
I figured providing the ungolfed version would help people understand the algorithm.  I have made marginal improvements since this version, but this should give you a jump start in understanding it:
import numpy as np

with open(&quot;test.txt&quot;,&quot;w&quot;) as fo:
    input = lambda: [[3,4], 2, 1]

    def bruteForce(array, size, number):
        for i in range(array.shape[0]-size+1):
            for j in range(array.shape[1]-size+1):
                if np.all(array[i:i+size,j:j+size] == 0):
                    temp = np.copy(array)
                    temp[i:i+size,j:j+size] = number
                    yield temp, (i+(size+1)/2,j+(size+1)/2)
                    
    #y, x because of how numpy works
    #euclid = https://codegolf.stackexchange.com/a/70536/67929
    def recursiveCheck(array,presents,presentLocations=[],number=1,currentMin=[10000,[]]):
        global rest, size
        if (len(presents) == 0):
            centerI = sum(a*b[0] for a,b in zip(rest,presentLocations)) / sum(rest)
            centerJ = sum(a*b[1] for a,b in zip(rest,presentLocations)) / sum(rest)
            euclid = lambda a,b:sum((d-e)**2for d,e in zip(a,b))**.5
            centerDelta = euclid([centerI,centerJ],[(a+1)/2 for a in array.shape])
            if centerDelta &lt;= currentMin[0]:
                currentMin[0] = centerDelta
                currentMin[1] = array
                fo.write(repr(array) + &quot;\n&quot;)
                fo.write(repr(presentLocations) + &quot;\n&quot;)
                fo.write(repr([centerI,centerJ]) + &quot;\n&quot;)
                fo.write(repr(centerDelta) + &quot;\n&quot;)
            return
        for potentialPresentArray, newPresentLocation in bruteForce(array,presents[0],number):
            recursiveCheck(potentialPresentArray,presents[1:],presentLocations+[newPresentLocation],number+1,currentMin)

    suMatrix = input()
    size = suMatrix[0]
    rest = suMatrix[1:]
    recursiveCheck(np.zeros(size),rest)


Credits
I used a response by CarpetPython to aide me in creating a euclidean distance calculation in a minimal amount of bytes.
Mr. Xcoder reduced it by 28 bytes.