5
$\begingroup$

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);
(*Earth ejection angle, in radians*)
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},
      {t, p}, Method -> "NelderMead"
     ]

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}}.

enter image description here

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).

Improve this question
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – MarcoB
    Commented Aug 9, 2015 at 1:45
  • $\begingroup$ 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. $\endgroup$
    – indigoblue
    Commented Aug 9, 2015 at 2:30
  • $\begingroup$ 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 $\endgroup$
    – indigoblue
    Commented Aug 9, 2015 at 2:43
  • $\begingroup$ 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? $\endgroup$
    – MarcoB
    Commented Aug 10, 2015 at 16:41
  • $\begingroup$ 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"] $\endgroup$
    – indigoblue
    Commented Aug 10, 2015 at 20:53

2 Answers 2

Reset to default
2
+50
$\begingroup$

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}, 
       ContourShading -> {White, Green}], 
     Graphics[{PointSize[Large], Red, Point[{rm[[2, 1, 2]], rm[[2, 2, 2]]}], Blue, 
       Point[{rmc[[2, 1, 2]], rmc[[2, 2, 2]]}]}]]

enter image description here

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.

Improve this answer
$\endgroup$
14
  • $\begingroup$ 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? $\endgroup$
    – indigoblue
    Commented Aug 16, 2015 at 16:57
  • $\begingroup$ @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.. $\endgroup$
    – bbgodfrey
    Commented Aug 16, 2015 at 17:31
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Commented Aug 16, 2015 at 18:08
  • $\begingroup$ 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. $\endgroup$
    – indigoblue
    Commented Aug 17, 2015 at 20:10
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Commented Aug 17, 2015 at 22:08
3
$\begingroup$

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}, 
             Method -> "NelderMead"]

(* {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]]]}]

Mathematica graphics

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]

Mathematica graphics

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.