ParametricNDSolve[ ] and NMinimize[ ] in Combination

Question

I have the following code that is used to simulate a spacecraft orbiting Mars, and uses ParametricNDSolve to vary the spacecraft's initial conditions using a parameter value p:

G = 6.672*10^-11;
m[1] = 6.4185*10^23;
m[2] = 1;
p[1] = {1000000, 1000000};
p[2] = {0, 0};
v[1] = {0, 0};
v[2] = {0, 2500};
tmax = 1000;

soln = ParametricNDSolve[{
    xm''[t] == -(G m[2] (xm[t] - xsc[t]))/((xm[t] - xsc[t])^2 + (ym[t] - ysc[t])^2)^(3/2),
    ym''[t] == -(G m[2] (ym[t] - ysc[t]))/((xm[t] - xsc[t])^2 + (ym[t] - ysc[t])^2)^(3/2),
    xsc''[t] == -(G m[1] (xsc[t] - xm[t]))/((xsc[t] - xm[t])^2 + (ysc[t] - ym[t])^2)^(3/2),
    ysc''[t] == -(G m[1] (ysc[t] - ym[t]))/((xsc[t] - xm[t])^2 + (ysc[t] - ym[t])^2)^(3/2),

    xm[0] == p[1][[1]], ym[0] == p[1][[2]], xsc[0] == p[2][[1]], 
    ysc[0] == p[2][[2]], xm'[0] == v[1][[1]], ym'[0] == v[1][[2]], 
    xsc'[0] == v[2][[1]] + 1000 p, ysc'[0] == v[2][[2]] + 1000 p},

   {xm[t], ym[t], xsc[t], ysc[t]}, {t, 0, tmax}, {p}, 
   Method -> "StiffnessSwitching", AccuracyGoal -> 18, 
   PrecisionGoal -> 18, MaxSteps -> 10000000];

   ParametricEscapeTrajectories[t_] = ParametricPlot[Evaluate[Table[{xsc[t][p], ysc[t][p]} /. soln, {p, 0, 5,1}]], {t, 0, tmax}, Prolog -> {Red, Disk[{p[1][[1]], p[1][[2]]}, 50000]}]

Below is an ensemble of results created by varying the parameter p from 0 to 5 as shown by the ParametricEscapeTrajectories function.

What I'm hoping is to be able to do something similar to what was done at this link, where the NMinimize function was used to find the closest approach distance between Mars and the spacecraft (the so-called periapse radius). However, this time I'm hoping to be able to find the value of closest approach, time of closest approach AND the value of the parameter p that gives me the closest approach out of all possible parameter values between 0 and 5. So, for example, it may look something like this:

{71729.9, {t -> 118.096}, {p -> 5}}

So far I've tried things such as the following, but unfortunately have not been successful:

NMinimize[Join[Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] /. soln, {t, 0, tmax}, {p, 0, 5, 1}], {t, p}]

Any help would be greatly appreciated, thanks very much.

Answer 1

NMinimize[{Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] /.  soln, 
           0 < t < tmax && 0 < p < 5}, 
           {t, p}, Method -> "SimulatedAnnealing"]

(* {71729.9, {t -> 118.095, p -> 5.}} *)

Edit

Please note that your system is more elegant and easier to understand when written in vector form:

G = 6.672*10^-11;
m[1] = 6.4185*10^23;
m[2] = 1;
p[1] = {1000000, 1000000};
p[2] = {0, 0};
v[1] = {0, 0};
v[2] = {0, 2500};
tmax = 1000;

soln = ParametricNDSolveValue[{
   posM''[t]  == -G m[2] (posM[t] - posSc[t])/  Tr[(posM[t] - posSc[t])^2]^(3/2),
   posSc''[t] == - m[1]/m[2] posM''[t],
   posM[0]  == p[1],  posSc[0] == p[2],
   posM'[0] == v[1],  posSc'[0] == v[2] + 1000 r },
   {posM, posSc}, {t, 0, tmax}, {r}, 
   Method -> "StiffnessSwitching", AccuracyGoal -> 18, 
   PrecisionGoal -> 18, MaxSteps -> 10000000]

ParametricPlot[Evaluate[Table[Last@Through[soln[p][t]], {p, 0, 5, 1}]], {t, 0, tmax}]