# NMinimize Constraints on Interpolating Functions

I have the following piece of code and am trying to use NMinimize (the last line in the following chunk of code) to minimize Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] with the constraint that this value must be greater than, say, 10000000. That is, Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}]>10000000:

Remove["Global*"]
(*Gravitational Constant*)
G = 6.672*10^-11;
(*Simulation running time, in seconds (1 day = 86400 seconds)*)
tmax = (1000) (86400);
(*Spacecraft time-of-flight between Earth and Mars, in seconds (1 day = 86400 seconds)*)
TOF = (254) (86400);
(*Mass of Sun, Earth, Mars (from JPL's Horizons ephemeris), in kilograms*)
ms = 1.988544*10^30 ;(*Mass of Sun*)
me = 5.97219*10^24 ;(*Mass of Earth*)
mm = 6.4185*10^23 ;(*Mass of Mars*)
(*Planetary radii of Sun, Earth and Mars (from JPL's Horizons ephemeris), in metres*)
r[0] = 6.963*10^8 ;(*Mean radius of Sun*)
r[1] = 6.37101 *10^6;(*Mean radius of Earth*)
r[2] = 3.3899*10^6 ;(*Mean radius of Mars*)
(*Heliocentric positions of Earth and Mars (from JPL's Horizons ephemeris) on 26 November 2011, in metres*)
pe = 149597870700 {4.416639858432274*10^-1, 8.798967504313304*10^-1} ;
pm = 149597870700 {-8.159382724017646*10^-1, 1.414986880765974*10^0};
(*Heliocentric velocities of Earth and Mars (from JPL's Horizons ephemeris) on 26 November 2011, in metres*)
ve = 149597870700/86400 {-1.563974110293042*10^-2, 7.690252775639107*10^-3};
vm = 149597870700/86400 {-1.160326991502370*10^-2, -5.778933879736245*10^-3} ;
(*Hyperbolic excess velocity, in metres per second. Difference between spacecraft's Heliocentric velocity vLambert[1] and Earth's Heliocentric velocity ve*)
vinf = Norm[{-29134.32, 15779.21} - ve];
(*Earth orbital radius of departure hyperbola, in metres*)
rp = (r[1] + 300000);
e = 6.37526;

(*Interplanetary transfer trajectory model*)
(*Models the motion of Earth (xe, ye), Mars (xm, ym) and spacecraft (xsc, ysc) in a Cartesian coordinate system*)
Soln = ParametricNDSolve[{
xe''[t] == -((G ms xe[t])/(xe[t]^2 + ye[t]^2)^(3/2)),
ye''[t] == -((G ms ye[t])/(xe[t]^2 + ye[t]^2)^(3/2)),
xm''[t] == -((G ms xm[t])/(xm[t]^2 + ym[t]^2)^(3/2)),
ym''[t] == -((G ms ym[t])/(xm[t]^2 + ym[t]^2)^(3/2)),
xsc''[t] == -((G ms xsc[t])/(xsc[t]^2 + ysc[t]^2)^(3/2)) - (G me (xsc[t] - xe[t]))/((xsc[t] - xe[t])^2 + (ysc[t] - ye[t])^2)^(3/2) - (G mm (xsc[t] - xm[t]))/((xsc[t] - xm[t])^2 + (ysc[t] - ym[t])^2)^(3/2),
ysc''[t] == -((G ms ysc[t])/(xsc[t]^2 + ysc[t]^2)^(3/2)) - (G me (ysc[t] - ye[t]))/((xsc[t] - xe[t])^2 + (ysc[t] - ye[t])^2)^(3/2) - (G mm (ysc[t] - ym[t]))/((xsc[t] - xm[t])^2 + (ysc[t] - ym[t])^2)^(3/2),

xe[0] == pe[[1]], ye[0] == pe[[2]], xm[0] == pm[[1]],
ym[0] == pm[[2]], xsc[0] == pe[[1]] + (r[1] + 300000) Cos[e p],
ysc[0] == pe[[2]] + (r[1] + 300000) Sin[e p], xe'[0] == ve[[1]],
ye'[0] == ve[[2]], xm'[0] == vm[[1]], ym'[0] == vm[[2]],
xsc'[0] == ve[[1]] - Sqrt[vinf^2 + (2 G me)/rp] Sin[e p],
ysc'[0] == ve[[2]] + Sqrt[vinf^2 + (2 G me)/rp] Cos[e p]}, {xe[t],
ye[t], xm[t], ym[t], xsc[t], ysc[t]}, {t, 0, tmax}, {p},
AccuracyGoal -> 8, PrecisionGoal -> 8,
Method -> "StiffnessSwitching", MaxSteps -> 10000000]

(*Finding value for parameter p that will fix the Earth ejection angle and produce the desired orbital radius upon Mars intercept*)
(*Constrained spacecraft to target an orbit of 300km above Mars*)
(*Constrained time search to be 10 days before and 10 days after "ideal" intercept time TOF*)
rm = NMinimize[{Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] /. Soln, TOF - 10 (86400) < t < TOF + 10 (86400) && 0.8 < p < 1}, {t, p}, Method -> "NelderMead"]


I thought that using may have worked:

rm = NMinimize[{
Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] /. Soln,
Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] > 10000000 /. Soln,
TOF - 10 (86400) < t < TOF + 10 (86400) && 0.8 < p < 1},
]


However, I seem to be getting syntax errors when adding in

Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] > 10000000 /. Soln

which I think is due to the fact that I'm not using /.Soln correctly.

Would anyone know how I can apply the desired constraint without getting such errors? Thanks very much.

EDIT:

Since the above bit of code it quite unwieldy, I've added a simpler example showing what I'd like to do. Suppose we have the following numerical output of a differential equation in the form of an interpolating function:

DE = Flatten[NDSolve[{u''[t] + u[t] == 0, u[0] == 0, u'[0] == 1}, u, {t, 0, 5}]]
Plot[u[t] /. DE, {t, 0, 5}]


And lets say I want to minimize the function u[t] on the interval 0<t<3, but also want u[t]>0.5. Then, using /.DE to apply the values of the interpolating function DE to u[t] I'd have

NMinimize[{u[t] /. DE, u[t] > 0.5 /. DE, 0 < t < 3}, t]


which gives me {0.5, {t -> 2.61799}}.

I'm trying to do something very similar to this, in the original code I posted. Except instead of minimizing u[t], I'm minimizing Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] and instead of requiring that u[t]>0.5, I would like Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] > 1000000. The only problem is that, following the same syntax procedures as I did on the simple example, I seem to be getting errors on the original code when trying to apply the additional constraint (the code gives no errors, however, without the additional constraint applied).

• Your code is too large and unwieldy. Please do some troubleshooting yourself first, to try and isolate the cause of your trouble, and then post a minimal example reproducing your problem. – MarcoB Aug 9 '15 at 1:45
• The issue is in the very last line of the code and everything above it can essentially be ignored. I put the rest of the code in so that someone running it would get results out. – indigoblue Aug 9 '15 at 2:30
• The cause of the trouble is in NMinimize. I'd like to put an additional constraint where Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}]>10000000, but cannot figure out how to do this since I get syntax errors. Every other piece of the code is working as planned – indigoblue Aug 9 '15 at 2:43
• So I ran your complex code chunk, and it runs without errors; the output is {1768.23, {t -> 2.10981*10^7, p -> 0.938622}}. I didn't get any syntax errors etc. Am I missing something? – MarcoB Aug 10 '15 at 16:41
• Thanks for taking a look MarcoB. Yup, the code currently runs without any errors, but I'm trying to add an additional constraint to NMinimize where Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}]>10000000 using the following rm = NMinimize[{Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] /. Soln, Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}]>10000000/.Soln, TOF - 10 (86400) < t < TOF + 10 (86400) && 0.8 < p < 1}, {t, p}, Method -> "NelderMead"] – indigoblue Aug 10 '15 at 20:53

As noted by the OP, adding the constraint,

Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] > 10000000 /. Soln


triggers error messages, probably because p typically must be given a numerical value before t is given one in expressions such as xm[t][p]. See the documentation for ParametricFunction. This can be accomplished by defining

g[t0_?NumericQ, p_?NumericQ] := Module[{tem =
Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] /. Soln}, tem /. t -> t0]


in which case the constrained expression, renamed rmc, becomes

rmc = NMinimize[{g[t, p], g[t, p] > 10000000, TOF - 10 (86400) < t < TOF + 10 (86400)
&& 0.8 < p < 1}, {t, p}, Method -> "NelderMead"]
(* {1.*10^7, {t -> 2.10972*10^7, p -> 0.938692}} *)


although with the inconsequential warning,

NMinimize::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations. >>


More insight can be obtained by plotting the curve in t - p space for which g[t, p] == 10000000, along with the points defined by rm and rmc.

Show[ContourPlot[g[t, p], {t, TOF - 9.85 (86400), TOF - 9.75 (86400)},
{p, .9385, .9387}, FrameLabel -> {t, p}, Contours -> {10000000},
Graphics[{PointSize[Large], Red, Point[{rm[[2, 1, 2]], rm[[2, 2, 2]]}], Blue,
Point[{rmc[[2, 1, 2]], rmc[[2, 2, 2]]}]}]]


The White region corresponds to Norm[...] < 10000000, which is excluded by the added constraint. The value of rmc (Blue) corresponds to but one point on a closed curve, Norm[...] == 10000000, and the value of rm (Red) to the true minimum computed in the question.

• Thanks bbgodfrey, I see that this did not work as I had hoped. I thought that, if I wanted to find the value of p that would minimize the distance between the spacecraft and Mars, but also added an additional constraint so that the spacecraft would not go below a certain value (in this case 10000000), NMinimize would find the global minimum to be the spacecraft's periapse radius relative to Mars and would not let this periapse radius go below 10000000. Would the "minimum of all minimums" for the above set of minimums be the periapse radius? – indigoblue Aug 16 '15 at 16:57
• @indigoblue To obtain what you describe in your last comment, use rm = NMinimize[{Norm[{xm[t][p], ym[t][p]} + {ysc'[t][p], -xsc'[t][p]} (r[2] + 10000000)/Norm[{xsc'[t][p], ysc'[t][p]}] - {xsc[t][p], ysc[t][p]}] /. Soln, TOF - 10 (86400) < t < TOF + 10 (86400) && 8/10 < p < 1}, {t, p}, Method -> "NelderMead"]. The calculation described in the question above instead yields all values of p for which the trajectory passes at 10000000 or closer to Mars, and the corresponding ts at which the distance is10000000, one while approaching and the other while departing.. – bbgodfrey Aug 16 '15 at 17:31
• @indigoblue Nonetheless, the plot in my answer above does give the values of p and t for which the closest approach of the trajectory is 10000000. They are the maximum and minimum p values on the curve, and the corresponding t values. One corresponds to passing Mars on the inside of its orbit, and the other on the outside. – bbgodfrey Aug 16 '15 at 18:08
• Do you possibly have a link to the mathematics used for the {ysc'[t][p], -xsc'[t][p]} (r[2] + 10000000)/Norm[{xsc'[t][p], ysc'[t][p]}] part? I see you use a unit vector formulation, but don't fully understand why {ysc'[t][p], -xsc'[t][p]} makes the spacecraft trajectory normal to a radial vector emanating from Mars' centre. – indigoblue Aug 17 '15 at 20:10
• @indigoblue I do not have a reference. Instead, I derived the expression by noting that a vector from the center of Mars to the trajectory must be perpendicular to the trajectory at the point of closest approach. See also the discussion in my answer to your earlier question. – bbgodfrey Aug 17 '15 at 22:08

As the Norm is positive, you don't really need the constraint to determine one minimum for p and t:

nm = NMinimize[{
(Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] - 10^7)^2 /. Soln,
TOF - 10 (86400) < t < TOF + 10 (86400) && 0.8 < p < 1}, {t, p},

(* {0.000155262, {t -> 2.10956*10^7, p -> 0.938678}} *)


Testing;

ContourPlot[ Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] /. Soln,
{t, TOF - 9.85 (86400), TOF - 9.75 (86400)}, {p, .9385, .9387},
Contours -> {10^7},
Epilog -> {Red, PointSize[Large], Point[{t, p} /. nm[[2]]]}]


You may find the parametric interpolation for the whole minimal contour:

cp = ContourPlot[
Norm[{xm[t][p], ym[t][p]} - {xsc[t][p], ysc[t][p]}] /. Soln,
{t, TOF - 9.85 (86400), TOF - 9.75 (86400)}, {p, .9385, .9387},
Contours -> {10^7}, ContourShading -> False, Frame -> False];

f = Interpolation /@ Transpose@First[Cases[cp // Normal, Line[x__] :> x, Infinity]]

ParametricPlot[Through[f[t]], {t, 1, 327}, AspectRatio -> 1]
`