Degenerate Conic

Minimization Without Derivatives

The FMIN function at Netlib is based on Richard Brent's classic localmin algorithm from [1]. The method uses a combination of golden section search and successive parabolic interpolation to compute the minimum of a 1D function without the necessity of evaluating any derivatives. There are Fortran 66 and ALGOL 60 versions of the method in the reference. I also have a modernized version here. Let's show an example of how to use it for a simple orbital mechanics application.

Problem setup

Let's solve a trivial problem in orbital mechanics: computing the time of closest approach to a planet for a conic orbit. This is known as periapsis (or, for the Earth, perigee). We don't really need to use a minimizer to compute this, but it's a problem with a known solution that we can test the minimizer on.

We need various algorithms to do this:

• A Kepler propagator and various utility routine such as Gooding's universal element transformations (we will use the Fortran Astrodynamics Toolkit)
• A minimizer (we will use the modernized fmin)

Given an initial elliptical orbit state, the Gooding routine pv3els will compute the time since last periapsis passage, so we have the "true" value to compare the results against. The function we pass to fmin to minimize will simply propagate the state (the independent variable is \(\Delta t\)) and return the radius magnitude (\(r\)). The point where the radius magnitude is minimized is the periapsis point we are looking for. The orbital period is used to set the upper time bound for the minimizer (we know there is only one periapsis passage per orbital period).

Example code

In our Fortran Package Manager manifest file, we define the two external dependencies we need:

[dependencies]
fmin = { git="https://github.com/jacobwilliams/fmin.git" }
fortran-astrodynamics-toolkit = { git="https://github.com/jacobwilliams/fortran-astrodynamics-toolkit.git" }

And the main program looks like this:

program main
  use fmin_module, wp => fmin_rk
  use fortran_astrodynamics_toolkit

  implicit none

  real(wp) :: p !! semiparameter [km]
  real(wp) :: period !! orbital period [sec]
  real(wp) :: dt !! time from initial state to periapsis [sec]
  real(wp),dimension(3) :: r, v !! position and velocity vectors [km] and [km/s]
  real(wp),dimension(6) :: rv0 !! initial state vector [km, km/s]
  real(wp),dimension(6) :: e0 !! Gooding orbital elements of the initial state
  integer :: n_evals = 0 !! number of function evaluations

  ! initial conditions (orbital elements):
  real(wp),parameter :: mu   = 398600.4418_wp  !! Earth grav. parameter [km^3/s^2]
  real(wp),parameter :: a    = 7000.0_wp       !! semi-major axis [km]
  real(wp),parameter :: ecc  = 0.2_wp          !! eccentricity
  real(wp),parameter :: inc  = 30.0_wp         !! inclination [deg]
  real(wp),parameter :: raan = 40.0_wp         !! right ascension of ascending node [deg]
  real(wp),parameter :: aop  = 60.0_wp         !! argument of perigee [deg]
  real(wp),parameter :: tru  = -45.0_wp        !! true anomaly [deg]
  real(wp),parameter :: tol  = 1.0e-6_wp       !! tolerance for fmin

  ! get the initial state vector and time from periapsis for the initial state:
  p = a * (1.0_wp - ecc**2)
  period = orbit_period(mu,a)
  call orbital_elements_to_rv(mu, p, ecc, inc, raan, aop, tru, r, v)
  rv0 = [r, v]
  call pv3els (mu, rv0, e0) ! e0(6) is the time from periapsis of the initial state

  ! call the minimizer:
  dt = fmin(func,0.0_wp,period,tol)

  ! print results
  write(*,'(A,F12.6,A)') 'DT to periapsis from fmin        = ', dt, ' sec'
  write(*,'(A,F12.6,A)') 'True initial time from periapsis = ', e0(6), ' sec'
  write(*,'(A,E12.3,A)') 'Error                            = ', dt + e0(6), ' sec'
  write(*,'(A,I0,A)')    'number of function evals         = ', n_evals, ' '

  contains

  function func(x) result(f)
    !! the function is to propagate from the initial state
    !! by the dt and compute the radius value
    real(wp),intent(in) :: x  !! indep. variable (dt)
    real(wp)            :: f  !! function value `f(x)` - radius magnitude
    real(wp),dimension(6) :: rvf
    call propagate(mu, rv0, x, rvf)
    f = norm2(rvf(1:3))
    n_evals = n_evals + 1
  end function func

end program main

Easy peasy!

Results

The result of running this program is:

DT to periapsis from fmin        =   652.983248 sec
True initial time from periapsis =  -652.983240 sec
Error                            =    0.875E-05 sec
number of function evals         = 12

So, it works! Using 12 function evaluations (i.e., kepler propagations of the orbit), the method has located the periapsis time to within about 8 \(\mu s\). This sort of approach can also be used for less-trivial problems where we don't have a good analytical solution, and/or derivatives are unavailable or hard to obtain.

References

1. R. Brent, "Algorithms for Minimization Without Derivatives", Prentice-Hall, 1973.

Sign, Sign, Everywhere a Sign

Reference [1] describes a numerical method for computing the derivative of a 1D function, using Neville's algorithm to extrapolate from a sequence of simple polynomial approximations based on interpolating points within specific bounds of a given point. The original code in the published paper in 1980 was written in ALGOL 60. It seems as if a Fortran 77 translation was produced by David Kahaner at NIST circa 1989. I found the file on this NIST server where it has been sitting quietly since 1992 (originally this was an ftp server when I first found it). It's an interesting algorithm, and I included a modernized version in my NumDiff numerical differentiation library. I also note that a version of the code is embedded within the Dataplot package.

But, there is a problem!

Here is a line of code in the original faccuracy function:

IF 16*ABS(H1)>ABS(H0) THEN H1:=SIGN(H1)*ABS(H0)/16;

Now, I know nothing about ALGOL 60, except that apparently John Backus (the creator of Fortran) was one of the committee members that developed it. Notice the SIGN function in this code. This line was translated into Fortran 77 (where the function is renamed as FACCUR) as:

IF(16.*ABS(H1) .GT. ABS(H0)) H1 = SIGN(H1,1.)*ABS(H0)/16.

They are very similar, but note the SIGN function again. In Fortran, the SIGN function (added in Fortran 77) is a little different since it has two arguments. While more flexible, it is a never-ending source of confusion, and I always have to look up the meaning of the two arguments. The following table summarizes the difference between the two:

So, do you see the bug? SIGN(H1,1.) doesn't return the same thing that SIGN(H1) did in the original code. SIGN(H1,1.) returns the value of H1 with the sign of 1.0, which is always +H1, not the original intent at all! It really should have been translated to SIGN(1.0,H1), to give the +1/-1 of the original code. So, this bug is over 30 years old. Technically, the 0 case is still not handled the same, but that's a degenerate case that would never happen or give meaningful results in this context. But to be totally faithful we would need to use a function like this:

elemental real(wp) function sign_algol(x)
real(wp),intent(in) :: x
sign_algol = merge(0.0_wp, sign(1.0_wp,x), x==0.0_wp)
end function sign_algol

Other Languages

Python doesn't have a built-in sign function, but the numpy one behaves exactly like the ALGOL one. The advantage of the Fortran one is that it does support signed zeros, which the Numpy one does not appear to support (for that you need to use copysign). C++ also has a copysign function that does the same thing.

References

1. J. Oliver, "Algorithm 017: An algorithm for numerical differentiation of a function of one real variable", Journal of Computational and Applied Mathematics 6 (2) (1980) 145-160. Fortran 77 code from NIST.
2. J. W. Backus, F. L. Bauer, J. Green, C. Katz, J. McCarthy, A. J. Perlis, H. Rutishauser, K. Samelson, B. Vauquois, J. H. Wegstein, A. van Wijngaarden, M. Woodger Editor: P. Naur, "Revised Report on the Algorithmic Language ALGOL 60", Communications of the ACM, Volume 6, Issue 1 Pages 1 - 17, 01 January 1963.
3. SIGN — Sign copying function [gcc.gnu.org]
4. Numerical Differentiation, degenerateconic.com, Dec 04, 2016