Downhill Simplex Optimization Analysis

downhillSimplexOptimizationAnalysis.py

# ------------------------------------------------------------------------
# The following Python code is implemented 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.
# Please see the Programming note for how to get started.
# ------------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt


# ------------------------------------------------------------------------
# INPUT
# ------------------------------------------------------------------------

# Objective function
def F(x):

    x1 = x[0]
    x2 = x[1]

    k1 = 100.0
    k2 = 90.0
    Fx1 = 20.0
    Fx2 = 40.0

    return k1 * ( np.sqrt( x1**2 + (x2+1)**2 ) - 1 )**2 + k2 * ( np.sqrt( x1**2 + (x2-1)**2 ) - 1 )**2 - ( Fx1 * x1 + Fx2 * x2 )


# Start point
startPoint = np.zeros(2)

# ------------------------------------------------------------------------
# SIMPLEX OPTIMIZATION ALGORITHM
# ------------------------------------------------------------------------

# Algorithm parameters
alpha = 1.0
gamma = 2.0
rho = -0.5
sigma = 0.5
tolerance = 0.00001
maxSteps = 100

# Initialize the plot
plt.ion()
plt.title('Downhill Simplex Optimization')

# Build a matrix of x-vectors
numDVs = len(startPoint)
xVectors = np.zeros((numDVs + 1, numDVs))

# Let the start point be the first x-vector
xVectors[0] = startPoint

# Build a simplex around the start point, creating each vertex by taking a step along each DV axis
for i in range(1, numDVs + 1):

    for j in range(numDVs):

        if j == i - 1:

            xVectors[i, j] = startPoint[j] + alpha

        else:

            xVectors[i, j] = startPoint[j]

# Calculate F at each vertex
Fvalues = np.zeros(numDVs + 1)
for i in range(numDVs + 1):
    Fvalues[i] = F(xVectors[i])

# Start the loop
loopCounter = 0
convergence = False
while not convergence and loopCounter < maxSteps:

    # Increment the loop counter
    loopCounter += 1

    # Update plot
    dv1Min = np.min([xVectors[0, 0], xVectors[1, 0], xVectors[2, 0], xVectors[0, 0]])
    dv1Max = np.max([xVectors[0, 0], xVectors[1, 0], xVectors[2, 0], xVectors[0, 0]])
    dv2Min = np.min([xVectors[0, 1], xVectors[1, 1], xVectors[2, 1], xVectors[0, 1]])
    dv2Max = np.max([xVectors[0, 1], xVectors[1, 1], xVectors[2, 1], xVectors[0, 1]])
    if loopCounter == 1:
        allTimeDv1Min = dv1Min
        allTimeDv1Max = dv1Max
        allTimeDv2Min = dv2Min
        allTimeDv2Max = dv2Max
    else:
        if dv1Min < allTimeDv1Min:
            allTimeDv1Min = dv1Min
        if dv1Max > allTimeDv1Max:
            allTimeDv1Max = dv1Max
        if dv2Min < allTimeDv2Min:
            allTimeDv2Min = dv2Min
        if dv2Max > allTimeDv2Max:
            allTimeDv2Max = dv2Max

    dv1Data = [xVectors[0, 0], xVectors[1, 0], xVectors[2, 0], xVectors[0, 0]]
    dv2Data = [xVectors[0, 1], xVectors[1, 1], xVectors[2, 1], xVectors[0, 1]]
    plt.clf()
    plt.plot(dv1Data, dv2Data, 'ko-', linewidth=1.0)
    plt.axis([allTimeDv1Min, allTimeDv1Max, allTimeDv2Min, allTimeDv2Max])
    plt.grid(True)
    plt.show()
    plt.pause(0.3)

    # Reset recorders
    reflectionBool = False
    expansionBool = False
    contractionBool = False
    shrinkingBool = False

    # Get the index of the worst and second-worst point
    worstValue = -1e+100
    bestValue = 1e+100
    for i in range(numDVs + 1):

        if Fvalues[i] > worstValue:
            worstValue = Fvalues[i]
            worstIndex = i

        if Fvalues[i] < bestValue:
            bestValue = Fvalues[i]
            bestIndex = i

    secondWorstValue = bestValue
    for i in range(numDVs + 1):

        if Fvalues[i] < worstValue and Fvalues[i] > secondWorstValue:
            secondWorstValue = Fvalues[i]
            secondWorstIndex = i

    # Calculate the centroid, considering all points except the worst
    centroid = np.zeros(numDVs)
    for i in range(numDVs + 1):

        if i != worstIndex:

            for j in range(numDVs):
                centroid[j] += xVectors[i, j]

    for j in range(numDVs):
        centroid[j] = centroid[j] / numDVs

    # Reflection
    reflection = np.zeros(numDVs)
    for i in range(numDVs):
        reflection[i] = centroid[i] + alpha * (centroid[i] - xVectors[worstIndex, i])

    F_reflection = F(reflection)

    # If the reflection point is better than the second worst, but not better than the best, accept it by replacing the worst point
    if F_reflection < secondWorstValue and F_reflection > bestValue:

        for i in range(numDVs):
            xVectors[worstIndex, i] = reflection[i]
            Fvalues[worstIndex] = F_reflection
            reflectionBool = True

        print('\n'"Accepted a reflection with F=", F_reflection)

    # If the reflection point is the best of all, then try further expansion
    elif F_reflection < bestValue:

        # Expansion
        expansion = np.zeros(numDVs)
        for i in range(numDVs):
            expansion[i] = centroid[i] + gamma * (reflection[i] - centroid[i])

        F_expansion = F(expansion)

        # If the expansion point is better than the reflection point, accept the expansion by replacing the worst point
        if F_expansion < F_reflection:

            for i in range(numDVs):
                xVectors[worstIndex, i] = expansion[i]
                Fvalues[worstIndex] = F_expansion
                expansionBool = True

            print('\n'"Accepted an expansion with F=", F_expansion)

        # If the expansion point isn't better than the reflection point, accept the reflection point by replacing the worst point
        else:

            for i in range(numDVs):
                xVectors[worstIndex, i] = reflection[i]
                Fvalues[worstIndex] = F_reflection
                reflectionBool = True

            print('\n'"Accepted a reflection with F=", F_reflection)

    # If the reflection point is neither best nor better than the second best, do contraction
    else:

        # Contraction
        contraction = np.zeros(numDVs)
        for i in range(numDVs):
            contraction[i] = centroid[i] + rho * (xVectors[worstIndex, i] - centroid[i])

        F_contraction = F(contraction)

        # If the contraction point is better than the worst point, accept it by replacing the worst point
        if F_contraction < worstValue:

            for i in range(numDVs):
                xVectors[worstIndex, i] = contraction[i]
                Fvalues[worstIndex] = F_contraction
                contractionBool = True

            print('\n'"Accepted a contraction with F=", F_contraction)

        # If not, then shrink
        else:

            # Shrinking
            shrinkingBool = True
            for i in range(numDVs + 1):
                if i != bestIndex:
                    for j in range(numDVs):
                        xVectors[i, j] = xVectors[bestIndex, j] + sigma * (xVectors[i, j] - xVectors[bestIndex, j])

                    Fvalues[i] = F(xVectors[i])

            print('\n'"Did a shrink")

    # Use max difference in F-values as convergence criterion for now
    if np.abs(bestValue - worstValue) < tolerance:

        convergence = True

        print('\n'"Converged with objective", Fvalues[0], "and design variables:")

        for j in range(numDVs):
            print(xVectors[0, j])

        print('\n'"Pausing a few seconds before closing the plot...")
        plt.pause(5)