goldenSectionLineSearch()

goldenSectionLineSearch.py

# ------------------------------------------------------------------------
# The following function is implemented in Python by Professor Terje Haukaas
# at the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without any form of warranty.
#
# The following notation applies:
# F(x) = objective function, or merit function, relevant in optimization, not in root-finding
# f(x) = dF/dx = derivative of the objective function, i.e., the function whose root is sought
# h(x) = df/dx = d^2F/dx^2 = second-order derivative, i.e., Hessian of the objective function
# ------------------------------------------------------------------------

def goldenSectionLineSearch(F, lowerBound, upperBound, maxIterations, tolerance, plot):

    # Initiate plot if requested (plot = delay if positive)
    if plot > 0:

        plt.clf()
        plt.ion()
        plt.title('Golden Section Search')
        plt.grid(True)
        plt.autoscale(True)
        plt.ylabel('Objective Function')
        plt.xlabel('Design Variable')
        xData = []
        FData = []

    # Calculate the golden ratio
    goldenRatio = (1.0 + np.sqrt(5.0)) / 2.0

    # Initialize the points  x1 ---------- x2 ----- x3 ---------- x4
    x1 = lowerBound
    x4 = upperBound
    range = x4 - x1
    bigPortionOfRange = range / goldenRatio
    x2 = x4 - bigPortionOfRange
    x3 = x1 + bigPortionOfRange

    # Evaluate the function at those points
    Fx1 = F(x1)
    Fx2 = F(x2)
    Fx3 = F(x3)
    Fx4 = F(x4)

    # Add points to the plot
    if plot > 0:

        xData.append(x1)
        xData.append(x2)
        xData.append(x3)
        xData.append(x4)

        FData.append(Fx1)
        FData.append(Fx2)
        FData.append(Fx3)
        FData.append(Fx4)

    # Start the loop
    convergence = False
    counter = 0
    while not convergence:

        # Increment counter
        counter += 1

        # Exit if we are at max number of iterations
        if counter > maxIterations:
            print("Golden section search reached maximum number of iterations without convergence")
            optimum = 0
            break

        # If the function value at x3 is greater than at x2 then the solution is in the range x1 ------------------- x2 ---------- x3
        # hence those points are now mapped to the new standard form                          x1 ---------- x2 ----- x3 ---------- x4
        if Fx3 > Fx2:

            # New range
            range = x3 - x1
            bigPortionOfRange = range / goldenRatio

            # Renaming previously calculated points
            x4 = x3
            Fx4 = Fx3

            x3 = x2
            Fx3 = Fx2

            # New point in that range
            x2 = x4 - bigPortionOfRange
            Fx2 = F(x2)

            # Add point to the plot
            if plot > 0:

                xData.append(x2)
                FData.append(Fx2)
                plt.plot(xData, FData, 'ko', linewidth=1.0)
                plt.show()
                plt.pause(plot)

        # Otherwise the solution is in the range                      x2 ---------- x3 ------------------- x4
        # hence those points are now mapped to the new standard form  x1 ---------- x2 ----- x3 ---------- x4
        else:

            # New range
            range = x4 - x2
            bigPortionOfRange = range / goldenRatio

            # Renaming previously calculated points
            x1 = x2
            Fx1 = Fx2

            x2 = x3
            Fx2 = Fx3

            # New point in that range
            x3 = x1 + bigPortionOfRange
            Fx3 = F(x3)

            # Add point to the plot
            if plot > 0:

                xData.append(x3)
                FData.append(Fx3)
                plt.plot(xData, FData, 'ko', linewidth=1.0)
                plt.show()
                plt.pause(plot)

        # Check convergence
        if np.abs(x3 - x2) < tolerance:

            optimum = (x2 + x3) / 2.0
            convergence = True

    # Output
    print('\n'"Golden section search done after", counter, "steps with solution", optimum)

    # Hold the plot for a few seconds before proceeding
    if plot > 0:
        print('\n'"Pausing a few seconds before closing the plot...")
        plt.pause(0.5)

    return optimum