Geek Answers Handbook

Why does scipy.optimize.minimize (default) report success without moving with Skyfield?

scipy.optimize.minimizeusing the default method, returns the original value as a result without any error messages or warnings. Although the Nelder-Mead method proposed by this answer solves the problem, I would like to understand:

Why does the default method return an incorrect answer without warning the starting point as an answer - and is there a way I can protect against a “wrong answer without warning” to avoid this behavior in this case?

Note that the function separationuses the python Skyfield package to generate minimized values ​​that are not guaranteed to be smooth, which may be the reason that Simplex is better here.

Results:

test result: [ 2.14159739 ] 'correct': 2.14159265359 initial: 0.0

default result: [ 10000. ] 'correct': 13054 initial: 10000

Nelder-Mead Result: [ 13053.81011963 ] 'correct': 13054 : 10000

FULL OUTPUT using DEFAULT METHOD:
   status: 0
  success: True
     njev: 1
     nfev: 3
 hess_inv: array([[1]])
      fun: 1694.98753895812
        x: array([ 10000.])
  message: 'Optimization terminated successfully.'
      jac: array([ 0.])
      nit: 0

FULL OUTPUT using Nelder-Mead METHOD:
  status: 0
    nfev: 63
 success: True
     fun: 3.2179306044608054
       x: array([ 13053.81011963])
 message: 'Optimization terminated successfully.'
     nit: 28

Here is the full script:

def g(x, a, b):
    return np.cos(a*x + b)

def separation(seconds, lat, lon):
    lat, lon, seconds = float(lat), float(lon), float(seconds) # necessary it seems
    place = earth.topos(lat, lon)
    jd = JulianDate(utc=(2016, 3, 9, 0, 0, seconds))
    mpos = place.at(jd).observe(moon).apparent().position.km
    spos = place.at(jd).observe(sun).apparent().position.km
    mlen = np.sqrt((mpos**2).sum())
    slen = np.sqrt((spos**2).sum())
    sepa = ((3600.*180./np.pi) *
            np.arccos(np.dot(mpos, spos)/(mlen*slen)))
    return sepa

from skyfield.api import load, now, JulianDate
import numpy as np
from scipy.optimize import minimize

data = load('de421.bsp')

sun   = data['sun']
earth = data['earth']
moon  = data['moon']

x_init = 0.0
out_g = minimize(g, x_init, args=(1, 1))
print "test result: ", out_g.x, "'correct': ", np.pi-1, "initial: ", x_init    # gives right answer

sec_init = 10000
out_s_def = minimize(separation, sec_init, args=(32.5, 215.1))
print "default result: ", out_s_def.x, "'correct': ", 13054, "initial: ", sec_init

sec_init = 10000
out_s_NM = minimize(separation, sec_init, args=(32.5, 215.1),
                 method = "Nelder-Mead")
print "Nelder-Mead result: ", out_s_NM.x, "'correct': ", 13054, "initial: ", sec_init

print ""
print "FULL OUTPUT using DEFAULT METHOD:"
print out_s_def
print ""
print "FULL OUTPUT using Nelder-Mead METHOD:"
print out_s_NM
+2
source share
3 answers

1)

- ( "" ). .

.

BFGS , ( , , , ).

, . , , , , . , , , .

, Nelder-Mead --- , , . , , BFGS.

2)

: . , .

, , , .

, . "" , .

/ , , --- , , .

+3

by @pv. , Skyfield "", , , , , .

- JulianDate, , 40 5E-10 . , , JPL . , , , . - "" (, ) . .

discrete time in skyfield

from skyfield.api import load, now, JulianDate
import numpy as np
import matplotlib.pyplot as plt

t  = 10000 + np.logspace(-10, 2, 25)        # logarithmic spacing
jd = JulianDate(utc=(2016, 3, 9, 0, 0, t))

dt  = t[1:] - t[:-1]
djd = jd.tt[1:] - jd.tt[:-1]

t  = 10000 + np.linspace(0, 0.001, 1001)        # linear spacing
jd = JulianDate(utc=(2016, 3, 9, 0, 0, t))

plt.figure()

plt.subplot(1,2,1)

plt.plot(dt, djd)
plt.xscale('log')
plt.yscale('log')

plt.subplot(1,2,2)

plt.plot(t, jd.tt-jd.tt[0])

plt.show()
+1

print , . separation(), , :

def separation(seconds, lat, lon):
    print seconds
    ...

, Nelder-Mead , 500 , :

[ 10000.]
[ 10500.]
[ 11000.]
[ 11500.]
[ 12500.]
...

, , 500- , . 500 , 500 , 500 . , Nelder-Mead, , , .

:

[ 10000.]
[ 10000.00000001]
[ 10000.]

. 1--8 , , , - .

​​, : (a) (b) , , , - , . - :

out_s_def = minimize(separation, sec_init, args=(32.5, 215.1),
                     tol=1e-3, options={'eps': 500})

- , , , - : , .

, , 64- , , . : " ", , , , , , !

, . :

  • float(), . print : , float, NumPy:

    (array([ 10000.]), 32.5, 215.1)

    : , Skyfield separation_from(), , :

    sepa = mpos.separation_from(spos)
return sepa.degrees

  • , Skyfield , 1.0.

- ( , , topos , ):

ts = load.timescale()

...

def separation(tt_jd, lat, lon):
    place = earth.topos(lat, lon)
    t = ts.tt(jd=tt_jd)
    mpos = place.at(t).observe(moon).apparent()
    spos = place.at(t).observe(sun).apparent()
    return mpos.separation_from(spos).degrees

...

sec_init = 10000.0
jd_init = ts.utc(2016, 3, 9, 0, 0, sec_init).tt
out_s_def = minimize(separation, jd_init, args=(32.5, 215.1))

, - , ? - , :

print ts.tt(jd=out_s_def.x).utc_jpl()

['A.D. 2016-Mar-09 03:37:33.8224 UT']

, Skyfield - , PyEphem , SciPy PyEphem C. , , - , : , .

Perhaps I should explore the possibility of creating Time objects from their two floating point objects, so that I can imagine a few more seconds. I think AstroPy did this, and this is traditionally in astronomy programming.

+1
source

More articles:


All Articles