Constraints in classical mechanics I am self-studying classical mechanics and I have a couple of questions about constraints. Goldstein in his book Classical Mechanics writes in the beginning of Section 1.3 that:

It is an overly simplified view to think that all problems in mechanics are reduced to solving the set of differential equations:
  $$m_i \ddot{\bf r}_i ~=~ {\bf F}_i^{(e)} + \sum_j {\bf F}_{ji} \tag{$\dagger$}$$ where ${\bf F}_i^{(e)}$ denotes the net external force on particle $i$ and ${\bf F}_{ji}$ denotes the force exerted by particle $j$ on $i$ because one may need to take into account the necessary constraints for the system.

Then, he says that the constraints introduce two types of difficulties in solving mechanical problems: 
$(\text{I})$ $r_i$ are no longer independent.
$(\text{II})$ constraint forces are not known in general.
My questions are:
$(1)$ IF it is possible to identify all the constraint forces, then all the problems would be reduced to solving $(\dagger)$, where ${\bf F}_i^{(e)}$ includes all the constraints forces. Wouldn't it? If not, is there a constraint that cannot be translated into a corresponding constraint force? 
$(2)$ Isn't $(\text{II})$ in fact the only difficulty with solving mechanical problems? It seems to me that $(\text{I})$ is not a "difficulty" because provided we can identify all the constraint forces the fact that $r_i$ are not independent would be incorporated into the constraint forces that would appear in the equations of motion. Isn't the fact that $r_i$ are not independent nothing but that the equations of motion are coupled ODEs? 
 A: Well, I think the book is not sufficiently clear. In Classical Mechanics without constraints, everything reduces to solve a system of differential equations of the form:
$$\frac{d^2 \vec{x}}{dt^2}= \vec{G}\left(t, \vec{x}(t), \frac{d \vec{x}}{dt}(t)\right)\tag{1}$$
with given initial conditions
$$\vec{x}(t_0) = \vec{x}_0\:,\quad \frac{d \vec{x}}{dt}(t) = \vec{v}_0\:.\tag{2}$$
where $\vec{x}=\vec{x}(t) \in \mathbb R^{3N}$ encompasses all the positions of the $N$ points of the system at time $t$,
$$\vec{x} = (\vec{x}_1, \ldots, \vec{x}_N)\:.$$
Above, $\vec{x}_i$ is the position vector of the $i$-th point in the rest 3-space of a reference frame. The masses $m_1,\ldots, m_N$ of the points have been embodied in the function $\vec{G}= \vec{G}(t, \vec{x}, \vec{v})$.
Since (1) is in normal form,  if the function $\vec{G}$ is known and is sufficiently regular, say $C^2$, (it would be enough jointly continuous in all variables and locally Lipschitz in $(\vec{x}, \vec{v})$), then the Cauchy problem (1)+(2) admits a unique (maximal) solution in a (maximal) interval including $t_0$.
Physically speaking the functional form of $\vec{G}$ is known when it is a superposition of (classically) fundamental forces, like the gravitational one, Lorentz force, etc...
However life is not so easy because there are practical and usual situations where we do not know the functional form of some of the forces superposed to produce $\vec{G}$. This is the case of geometrically constrained points, where, instead of the functional form of the force necessary to make the constraint satisfied, only the analytic equation of the constraint is provided. In this case (1) has to be replaced by
$$\frac{d^2 \vec{x}_i}{dt^2}= \frac{1}{m_i}\vec{F}_i\left(t, \vec{x}(t), \frac{d \vec{x}}{dt}(t)\right)+ \vec{\phi}_i(t)\tag{1'}\quad i=1,\ldots, N$$
where $\vec{\phi}_i(t)$ is the force on the $i$-th point due to the geometric constraint at time $t$, whose functional form is completely unknown and therefore it is a further unknown of the problem exactly as the functions $\vec{x}_i$. 
A closer microscopic scrutiny would reveal that $\phi_i$ is of electrical nature, but it is irrelevant here. They also can be mathematically emulated by means of a suitable functional forms, e.g., of conservative forces, taking the form of the constraints into account. However, barring particular cases, these approaches do not make easier the solution of the problem of the constrained motion.
(1)' has to be accompanied by a set of equations describing the constraints
$$f_j(t, \vec{x}_1, \ldots, \vec{x}_N)=0, \quad j=1,\ldots, c\tag{1''}$$
with $3N-c >0$ and the functions $f_j$ have to be functionally independent (I will not insist here on this point).
The set of requirements (1')+(1'')+(2) (where (2) must be compatible with the constraints) are generally not yet enough to determine a unique solutions 
$\vec{x}_i= \vec{x}_i(t)$ together with the values of the functions $\vec{\phi}_i= \vec{\phi}_i(t)$ along every motion of the system.
The remaining necessary information is given by a constitutive equation about the constraints. This is a non-geometrical relation involving the unknown forces $\vec{\phi}_i$. You should know several cases. The most important are two. The standard friction relation about components of the $\vec{\phi}_i$ involving friction coefficients $\mu$ and the D'Alembert's characterization of the ideal constraints. The latter has as the simplest case a frictionless constraint, but also includes the rigidity constraint and the integrable rolling constraint.
The problem consisting of (1')+(1'')+(2)+D'Alembert's characterization of the ideal constraints always admits  a unique solution 
$$\vec{x}_i= \vec{x}_i(t)\:, \quad \vec{\phi}_i=\vec{\phi}_i(t)\:, i=1,\ldots, N$$
provided all known involved functions are sufficiently regular. This is the main result of Lagrangian formulation of classical mechanics. As a byproduct, the equation of motion are written using just $3N-c$ coordinates (the $c$ constraints as well as the unknown forces $\vec{\phi}_i$ disappear from the Lagrange equations) which can be chosen with a very great arbitrariness as you probably already know.

Addendum 1.
Regarding the charcterization of D'Alembert constarints it is like this. $$\sum_{i=1}^N \vec{\phi}_i \cdot \delta \vec{x}_i =0 \tag{D}$$
where everything is evaluated at every fixed time $t$, on every permitted configuration $\vec{x}$ and each vector $$\delta \vec{x} = (\delta\vec{x}_1, \ldots, \delta\vec{x}_N)$$ is every tangent vector to the submanifold of $\mathbb R^{3N}$ including all permitted configurations at time $t$.
If we introduce local free coordinates $q^1,\ldots, q^n$ ($n= 3N-c$) on this manifold, we have 
$$ \delta \vec{x}_i := \sum_{k=1}^n\frac{\partial \vec{x}_i}{\partial q^k} \delta q^k \tag{E}$$
for every choice of the numbers $\delta q^k \in \mathbb R$.
$\delta \vec{x}$ is sometimes called virtual displacement.
(The notation is horrible) $\delta q^k$ are not small, they are arbitrary. I know that some books write "infinitesimal", but it does not mean anything. The right-hand side of (E) (with $i$ fixed) is a generic tangent vector to the manifold of the constraints, the $\delta q^k$ are the generic scalar components of that vector with respect to the basis made of the $n$ vectors $\frac{\partial \vec{x}_i}{\partial q^1}\:, \ldots,\frac{\partial \vec{x}_i}{\partial q^n}$.  These vectors are linearly independent because the constraints are assumed to be functionally independent.
(D) Says that the set of  reactive forces $\phi_i$ is always orthogonal to the manifold of constraints in a generalized sense. This is the geometric meaning of the ideal contraints.
It is possible to prove that (D) includes the case of a frictionless geometric constraints (also when the geometric curves and surfaces change their shape in time), the rigidity constraint, any mixing of these two cases. The case of rolling constraint is also included provided it is integrable, i.e., it can be re-stated into a pure geometric fashion as it happens for a disk rolling on a fixed path. But it is not the case for a sphere rolling on a plane.   

Addendum 2.
Regarding the existence and uniqueness of the solution of the arising equations of motion. First of all, (D) and (E) permit us to restate everything in terms of free local coordinates $q^1,\ldots, q^n$. The obtained  equations of motion are the known general Euler-Lagrange equations
$$\frac{d}{dt} \frac{\partial T}{\partial \dot{q}} - 
 \frac{\partial T}{\partial q} = Q_k\tag{EL1}$$
$$\frac{dq^k}{dt} = \dot{q}^k\tag{EL2}$$
where, for $k=1,\ldots, n$
$$Q_k(t,q, \dot{q}) := \sum_{i=1}^N \vec{F}_i \cdot \frac{\partial \vec{x}_i}{\partial q^k} $$
and 
$$T(t,q,\dot{q}) = \sum_{i=1}^N \frac{1}{2} m_i \left(\frac{d\vec{x}_i}{dt}\right)^2\:.$$
It is possible to prove that, under our hypotheses of ideal geometric constraints 
$$T(t,q,\dot{q}) = \sum_{k,h=1}^m a_{hk}(t,q) \dot{q}^h \dot{q}^k  + \sum_{k=1}^m b_{k}(t,q) \dot{q}^k + c(t,q)$$
where the matrix $[a_{hk}]$ is non-singular and positive. As a matter of fact
$$\begin{align}
a_{hk}(t,q) &= \sum_{i=1}^N \frac{1}{2} m_i \frac{\partial \vec{x}_i}{\partial q^h} \cdot  \frac{\partial \vec{x}_i}{\partial q^k} \\
b_{k}(t,q) &= \sum_{i=1}^N m_i \frac{\partial \vec{x}_i}{\partial q^k} \cdot  \frac{\partial \vec{x}_i}{\partial t} \\
c(t,q) &= \sum_{i=1}^N \frac{1}{2} m_i \frac{\partial \vec{x}_i}{\partial t} \cdot  \frac{\partial \vec{x}_i}{\partial t}
\end{align}$$
Using this fact, a complicated computation proves that (EL1)+(EL2) reduces to a system of the from
$$\frac{d^2 q^k}{dt^2} = G^k\left(t, q, \frac{d q}{dt}  \right) \quad k=1,\ldots$$
This is a system of differential equations of the second order written into its normal form (all highest order derivatives appear separated from the other variables).
If the right-hand side is sufficiently regular (as is the case if the constraints and the functional form of the known forces are smooth), the general theorem of existence and uniqueness of the solutions for initial conditions $q^k(t_0)$, $\frac{dq^k}{dt}|_{t_0}$ applies.   
A: *

*Constraints can in principle be emulated by introducing appropriate generalized stiff springs. Then the positions ${\bf r}_i$ remain independent. We choose the spring constants $k_j$ so big (but finite $<\infty$) that the constraints are satisfied to any precision we may want. (And this approach may actually be closer to how a realistic system behaves.) 

*However, it is much simpler to solve a system with idealized constraints (with corresponding possible infinite constraint forces to rigorously impose these constraints).  Idealized constraints mean that the positions ${\bf r}_i$ are no longer independent, as Goldstein writes. Of course it's no fun to have infinite constraint forces around. However, one can often argue that constraint forces do no virtual work. This leads to d'Alembert's principle, cf. e.g. this Phys.SE post and links therein.
A: 1) It can be extremely difficult to translate a constraint into a force. Think for example about the classic problem of a bead sliding along a rail: how would you write the force corresponding to such a constraint? It is much easier to impose that the motion of the bead con only take place along the rail and adopt a different set of coordinates instead of the usual Cartesian coordinates. This is why it is usually much easier to solve such problems using the Lagrangian rather than the Newtonian approach. 
2) Yes, the equations would "only" become coupled. However, this is not just some minor disadvantage. It is much more difficult to solve a system of coupled equations rather than a set of uncoupled equations. Just think about a set of $N$ uncoupled harmonic oscillators: it is trivial to write down their equation of motion. But what if we couple some of them together? I assure you that the equation of motion can then become a real mess.
