bisectionLineSearch()

bisectionLineSearch.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 bisectionLineSearch(F, lowerBound, upperBound, maxIterations, tolerance, plot):

    # Initialize plot
    if plot > 0:
        plt.clf()
        plt.ion()
        plt.title('Bisection Search')
        plt.grid(True)
        plt.autoscale(True)
        plt.ylabel('Derivative of Objective Function')
        plt.xlabel('Design Variable')
        xData = []
        fData = []

    # Evaluate f at the lowerbound (the lower bound to trial point is used as reference interval to test where the solution is)
    f_lower = f(lowerBound, F)
    f_lowerNeedsEvaluation = False

    # Check that the solution is indeed between the lower and upper bounds
    f_upper = f(upperBound, F)
    if not f_lower * f_upper < 0.0:
        print("Bisection search did not find a solution between the lower and upper bounds")

    # Add points to the plot
    if plot > 0:
        xData.append(lowerBound)
        xData.append(upperBound)
        fData.append(f_lower)
        fData.append(f_upper)
        plt.plot(xData, fData, 'ko', linewidth=1.0)
        plt.show()
        plt.pause(plot)

    # 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("Bisection search reached the maximum number of iterations without convergence")
            optimum = 0
            break

        # Calculate trialPoint, the midpoint between the lower and upper bounds
        trialPoint = (upperBound + lowerBound) / 2.0

        # Evaluate f at trialPoint
        f_trial = f(trialPoint, F)

        # Add point to the plot
        xData.append(trialPoint)
        fData.append(f_trial)
        plt.plot(xData, fData, 'ko', linewidth=1.0)
        plt.show()
        plt.pause(plot)

        # Evaluate f at the lower bound if needed, and if so, add the point to the plot
        if f_lowerNeedsEvaluation:

            f_lower = f(lowerBound, F)

            # Add point to the plot
            if plot > 0:
                xData.append(lowerBound)
                fData.append(f_lower)
                plt.plot(xData, fData, 'ko', linewidth=1.0)
                plt.show()
                plt.pause(plot)

        # Check if trialPoint should replace the lower or the upper bound
        if f_lower * f_trial < 0.0:
            upperBound = trialPoint
            f_lowerNeedsEvaluation = False
        else:
            lowerBound = trialPoint
            f_lowerNeedsEvaluation = True

        # Check convergence
        if np.abs(upperBound - lowerBound) < tolerance:

            optimum = trialPoint
            convergence = True

    # Output
    print('\n'"Bisection 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