Partial Answer/Too long for comment. If you are just working in the plane, then intuitively you know already that a length-minimizer exists (so long as you allow self-intersections). Such a solution will be a geodesic. This means you know (ahead of designing your functional) that a solution won't have any curvature - so one ought to focus on the length minimzation.

From here, it seems to me (I have not done any detailed calculation) that some sort of reasonable convex penalization on your curves for not going through a given point $q$ will result in the right critical points. For example, if you work in a class of piecewise $C^1$ curves $\gamma :(0,1) \to \mathbb{R}^2$, then consider

$F(\gamma) = \text{length}(\gamma) + \inf_{x \in (0,1)}|p-q|^2$.

You can then calculate the Euler-Lagrange equation for this functional. You do need to justify differentiating through the infimum, though.

The more interesting question is that of visiting multiple points in some order. I have not thought about this.

Partial Answer/Too long for comment. If you are just working in the plane, then intuitively you know already that a length-minimizer exists (so long as you allow self-intersections). Such a solution will be a geodesic. This means you know (ahead of designing your functional) that a solution won't have any curvature - so one ought to focus on the length minimzation.

From here, it seems to me (I have not done any detailed calculation) that some sort of reasonable convex penalization on your curves for not going through a given point $q$ will result in the right critical points. For example, if you work in a class of piecewise $C^1$ curves $\gamma :(0,1) \to \mathbb{R}^2$, then consider

$F(\gamma) = \text{length}(\gamma) + \inf_{x \in (0,1)}|p-q|^2$.

You can then calculate the Euler-Lagrange equation for this functional. You do need to justify differentiating through the infimum, though.

The more interesting question is that of visiting multiple points in some order. I have not thought about this.