Consider a system of $N$ material points described by position vectors ${\bf x}_1, \ldots, {\bf x}_N$ in a reference frame ${\cal R}$. These position vectors are not free to assume any configuration in the physical space, but they are constrained to satisfy some constraints which, possibly, may depend on time,
$$f_j({\bf x}_1, \ldots, {\bf x}_N, t) =0\quad j=1,\ldots, c < 3N\:.\tag{1}$$
If these functions are smooth and satisfy a condition of functional independence (I do not want to enter into the details), we can choose $n:= 3N-c$ abstract coordinates $q^1,\ldots, q^n$ which can be used the embody the constraints into the formalism. This result holds locally around every admitted configuration and around a given time $t$.
As a matter of fact, we can locally (in space and time) represent the position vectors ${\bf x}_1, \ldots, {\bf x}_N\quad i=1,2,\ldots, N$ as known functions of the said free coordinates.
$${\bf x}_i = {\bf x}_i(q^1,\ldots, q^n,t)\tag{2}$$
When $q^1,\ldots, q^n, t$ varies in their domain (an open set in $\mathbb{R}^{n+1}$), the vectors ${\bf x}_i(q^1,\ldots, q^n,t)$ automatically satisfy the constraints (1). The admissible configuration are therefore determined by the free coordinates $q^1,\ldots, q^n$ at each time $t$.
Now we pass to the notion of virtual displacement compatible with the set of constraints (1). It is defined by fixing $t$ and computing the differential of the functions ${\bf x}_i$ as functions of the remaining variables. The virtual displacement of the system at time $t_0$ around a permitted configuration determined by $q_0^1,\ldots, q_0^n$ is the set of $N$ vectors in the real space
$$\delta {\bf x}_i = \sum_{k=1}^n \left.\frac{\partial {\bf x}_i}{\partial q^k}\right|_{(q^1_0,\ldots, q_0^n,t_0)}\delta q^k\:, \quad i=1,\ldots, N\:.$$
Above the numbers $\delta q^k\in {\mathbb R}$ are arbitrary, not necessarily "infinitesimal" (which does mean anything!).
Example
Consider a point of position vector ${\bf x}$, constrained to live on a circle of radius $r=\sqrt{1+ct^2}$, where $c>0$ is a known constant.
This circle is centered on the origin and stays in the plane $z=0$.
Here we have just two constraint
$$f_1({\bf x}, t)=0, \quad f_2({\bf x}, t)=0$$
where
$$f_1({\bf x}, t) := x^2+ y^2+ z^2 - (1+ct^2)\:, \quad
f_1({\bf x}, t) := z $$
if ${\bf x}= x{\bf e}_x+ y{\bf e}_y+ z{\bf e}_z$
Locally we can use, for instance, the coordinate $q^1=x$ to describe a portion of circle and we have
$$x = q^1\:, \quad y = \sqrt{(1+ct^2) - x^2}\:, \quad z=0\:.$$
The relations (2) here read
$${\bf x}(q^1,t)= q^1 {\bf e}_1+ \sqrt{(1+ct^2) - (q^1)^2}{\bf e}_2$$
The virtual displacements at time $t_0$ are the vectors of the form
$$\delta {\bf x} = \delta q^1 {\bf e}_1+ \frac{\delta q^1}{\sqrt{(1+ct_0^2) - (q^1)^2}}{\bf e}_2$$
for every choice of $\delta q^1$.
The geometric meaning of $\delta {\bf x}$ should be evident: it is nothing but a vector (of arbitrary length) tangent to the circle at time $t_0$ emitted by a configuration determined by the value $q^1$.
REMARK It is worth stressing that varying $t$, the circle changes! Virtual displacement are defined at given time $t$.
The discussed example is actually quite general. The virtual displacement are always vectors which are tangent to the manifold of admitted configurations at given time $t_0$.
The prosecution of Lagrangian approach (once assumed the validity of the postulate of ideal constraints and introducing interactions for instance defined by a Lagrangian) consists of finding the evolution of the system not in terms of curves $${\bf x}_i = {\bf x}_i(t)\:, \quad i=1\,\ldots, N$$ in the physical space. The motion is described directly in terms of free coordinates, i.e., curves $$q^k=q^k(t)\:, \quad k=1,\ldots, n$$
Just at this level it makes sense to introduce the notation $\dot{q}^k = \frac{dq^k}{dt}$, because here we have a curve $t$-parametrized describing the evolution of the system. Before this step $q^1,\ldots, q^n$ and $t$ are independent variables.
A posterori, if we have a description of the motion of the system in terms of free coordinates, we also have the representation in the physical space by composing these curves with the universal (= independent of any possible motion) relations (2),
$${\bf x}_i(t) = {\bf x}_i(q^1(t), \ldots, q^n(t), t)\:, \quad i=1\,\ldots, N\tag{3}$$
It is finally interesting to compare virtual displacements with real displacements when we have a motion $$q^k=q^k(t)\:, \quad k=1,\ldots, n\:.$$
In the physical space, the velocities with respect to the reference frame ${\cal R}$ are given by taking the derivative with respect to $t$ of (3),
$${\bf v}_i(t) = \sum_{k=1}^n \frac{\partial {\bf x}_i}{\partial q^k}\frac{dq^k}{dt} + \frac{\partial {\bf x}_i}{\partial t}\:.$$
An approximate displacement ascribed to an interval of time $\Delta t$ is
$$\Delta {\bf x}_i = {\bf v}_i(t)\Delta t = \sum_{k=1}^n \frac{\partial {\bf x}_i}{\partial q^k}\frac{dq^k}{dt}\Delta t + \frac{\partial {\bf x}_i}{\partial t}\Delta t\:,$$
which can be rephrased to
$$\Delta {\bf x}_i = \sum_{k=1}^n \frac{\partial {\bf x}_i}{\partial q^k}\Delta q^k + \frac{\partial {\bf x}_i}{\partial t}\Delta t\:. \tag{4}$$
This identity has to be compared with the definition of virtual desplacement
$$\delta {\bf x}_i = \sum_{k=1}^n \frac{\partial {\bf x}_i}{\partial q^k}\delta q^k\:.$$
Even if we choose $\delta q^k = \Delta q^k$, the right-hand sides are different in view of the term $ \frac{\partial {\bf x}_i}{\partial t}\Delta t$ which accounts for a part of the displacement, in real motion, due to the fact that constraints may depend on $t$ explicitly, as in the example above.
