Finding least-norm solutions of a linear Diophantine equation I am looking for an algorithm solving the following problem.

Let $\Vert x\Vert_{*} := \sqrt{x^\mathrm{T} D \,x}$, where $D$ is a diagonal matrix with positive diagonal entries.
Given $a\in \mathbb{N}^n$, $b \in \mathbb{N}$, find
$$\arg \min_{x \in \mathbb{N}^n} {\vert a^{\mathrm{T}} x - b\vert}$$
with minimal $\Vert x \Vert_{*}$.

Motivation
I am creating a tools website for tabletop role-playing game (TTRPG) players. In these games, players fight together against monsters, and one player (called "GM") has to ensure these fights are of appropriate difficulty. This is done using "experience points" (XP). XP is used to estimate the strength of monsters. Each monster awards a set amount of XP (tougher monsters award more XP) and each difficulty has a target XP number. When the players defeat the monsters they gain these XP, becoming stronger as a result. If the GM selects a few monsters, how many of each do we need to ensure a good battle?
For example, the GM selects goblins, hobgoblins, and bugbears (who award 1,2, and 4 XP) to be his monsters. As threshold they select 12 XP. We would have $a^T=[1,2,4], b=12.$ $x$ describes the number of monsters picked. It is clear that there are many ways to select monsters, but most of them are undesirable. We don't want a fight against just 12 goblins, or only 3 bugbears. My idea is giving each monster a somewhat arbitrary weight: $D=diag(1,6,12)$. With this we obtain the optimal $x^T=[4,2,1]$ (4 goblins, 2 hobgoblins and 1 bugbear), a well balanced encounter.
If the GM selects creatures such that we can't meet the threshold, we want to get as close as possible to it.
Why don't I use a greedy algorithm? Because I thought it would be cool and interesting not to.
 A: You seem to be juggling two criteria (getting close to the threshold and minimizing the weighted norm of $x$). There are several ways you could approach this via an optimization model. One would be to minimize a weighted combination of $\Vert x \Vert_*^2$ and $\vert a^Tx - b\vert$. (I squared the norm so that the objective function is quadratic and convex.) Another would be to minimize the squared norm subject to an arbitrary limit on the absolute deviation. Alternatively, you could minimize the absolute deviation subject to a limit on the squared norm. All three approaches should lead to tractible problems.
A: First,
$$\begin{array}{ll} \underset{\mathrm x \in \Bbb N^3}{\text{minimize}} & \left| \, {\rm a}^\top {\rm x} - b \,  \right|\end{array}$$
If the minimum is zero, then solve the following least-norm problem
$$\begin{array}{ll} \underset{\mathrm x \in \Bbb N^3}{\text{minimize}} & \| \mathrm x \|_*^2\\ \text{subject to} & \mathrm a^\top \mathrm x = b\end{array}$$
Using Python (NumPy + CVXPY),
import cvxpy as cp
import numpy as np

# parameters
a = np.array([1, 2, 4])
b = 12
D = np.diag([1, 6, 12])

# optimization variables
x = cp.Variable(3, integer=True)


# create and solve linear optimization problem over the integers
objective = cp.Minimize(cp.abs(a.T @ x - b))
prob      = cp.Problem(objective, [ x >= np.zeros(3) ])
solution  = prob.solve()

print("—————————————————————————————")
print("IP ——————————————————————————")
print("Minimum  = ", prob.value)
print("Solution = ",    x.value)


# create and solve quadratic optimization problem over the integers
objective2  = cp.Minimize(cp.quad_form(x, D))
constraints = [ a.T @ x - b ==  0,
                x >= np.zeros(3)  ]
prob2       = cp.Problem(objective2, constraints)
solution2   = prob2.solve(solver=cp.ECOS_BB)

print("—————————————————————————————")
print("MIQP ————————————————————————")
print("Minimum  = ", prob2.value)
print("Solution = ",     x.value)

which outputs
—————————————————————————————
IP ——————————————————————————
Minimum  =  0.0
Solution =  [0. 0. 3.]
—————————————————————————————
MIQP ————————————————————————
Minimum  =  51.99999927150701
Solution =  [4. 2. 1.]

Unfortunately, according to the CVXPY documentation, there are correctness issues with this solver.
