3

I am new to python and in learning stages. I wanted to implement Particle Swarm Optimization(PSO) algorithm which I did by taking help from on-line materials and python tutorials. In PSO, a simple calculus problem is inferred i-e 100 * ((y - (x2))2) + ((1 - (x2))2). This problem is defined in a fitness function.

def fitness(x, y):
    return 100 * ((y - (x**2))**2) + ((1 - (x**2))**2)

Now, I want to replace this simple calculus problem by simple first order Ordinary Differential Equation(ODE) by without changing existing function parameters (x,y) and want to return the value of dy_dx,y0 and t for further process.

# Define a function which calculates the derivative
def dy_dx(y, x):
    return x - y

t = np.linspace(0,5,100)
y0 = 1.0  # the initial condition
ys = odeint(dy_dx, y0, t)`

In python odeint function is used for ODE which requires three essential parameters i-e func/model, y0( Initial condition on y (can be a vector) and t(A sequence of time points for which to solve for y) Example of odeint parameters. I don't want to change its parameters because it will be difficult for me to make changes in algorithm.

For simplicity I pasted the full code below and my question is open to anyone if wants to modify the code with further parameters in General Best, Personal Best and r[i].

    import numpy as np
    from scipy.integrate import odeint
    import random as rand
    from scipy.integrate import odeint
    from numpy import array
    import matplotlib.pyplot as plt

    def main():
    #Variables
    n = 40
    num_variables = 2

    a = np.empty((num_variables, n))
    v = np.empty((num_variables, n))
    Pbest = np.empty((num_variables, n))
    Gbest = np.empty((1, 2))
    r = np.empty((n))

    for i in range(0, num_variables):
        for j in range(0, n):
            Pbest[i][j] = rand.randint(-20, 20)
            a[i][j] = Pbest[i][j]
            v[i][j] = 0

    for i in range(0, n):
        r[i] = fitness(a[0][i], a[1][i])

    #Sort elements of Pbest
    Order(Pbest, r, n)

    Gbest[0][0] = Pbest[0][0]
    Gbest[0][1] = Pbest[1][0]

    generation = 0

    plt.ion()
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.grid(True)

    while(generation < 1000):
        for i in range(n):
            #Get Personal Best
            if(fitness(a[0][i], a[1][i]) < fitness(Pbest[0][i], Pbest[1][i])):
                Pbest[0][i] = a[0][i]
                Pbest[1][i] = a[1][i]
            #Get General Best
            if(fitness(Pbest[0][i], Pbest[1][i]) < fitness(Gbest[0][0], Gbest[0][1])):
                Gbest[0][0] = Pbest[0][i]
                Gbest[0][1] = Pbest[1][i]
            #Calculate Velocity
            Vector_Velocidad(n, a, Pbest, Gbest, v)
        generation = generation + 1
        print 'Generacion: ' + str(generation) + ' - - - Gbest: ' +str(Gbest)

        line1 = ax.plot(a[0], a[1], 'r+')
        line2 = ax.plot(Gbest[0][0], Gbest[0][1], 'g*')

        ax.set_xlim(-10, 10)
        ax.set_ylim(-10, 10)

        fig.canvas.draw()

        ax.clear()
        ax.grid(True)

    print 'Gbest: '
    print Gbest

def Vector_Velocidad(n, a, Pbest, Gbest, v):
    for i in range(n):
        #Velocity in X
        v[0][i] = 0.7 * v[0][i] + (Pbest[0][i] - a[0][i]) * rand.random() * 1.47 + (Gbest[0][0] - a[0][i]) * rand.random() * 1.47
        a[0][i] = a[0][i] + v[0][i]

        v[1][i] = 0.7 * v[1][i] + (Pbest[1][i] - a[1][i]) * rand.random() * 1.47 + (Gbest[0][1] - a[1][i]) * rand.random() * 1.47
        a[1][i] = a[1][i] + v[1][i]

def fitness(x, y):
        return 100 * ((y - (x**2))**2) + ((1 - (x**2))**2)

def Order(Pbest, r, n):

    for i in range(1, n):
        for j in range(0, n - 1):
            if r[j] > r[j + 1]:
                #Order the fitness
                tempRes = r[j]
                r[j] = r[j + 1]
                r[j + 1] = tempRes

                #Order las X, Y
                tempX = Pbest[0][j]
                Pbest[0][j] = Pbest[0][j + 1]
                Pbest[0][j + 1] = tempX

                tempY = Pbest[1][j]
                Pbest[1][j] = Pbest[1][j + 1]
                Pbest[1][j + 1] = tempY

if '__main__' == main():
    main()
Usman YousafZai
  • 1,088
  • 4
  • 18
  • 44
  • 1
    You should be able to replace your `Ordenamiento_Burbuja` with `argsort` and some indexing – Eric Jan 05 '18 at 13:59
  • @Eric ... I have tried it from every angle but I am not able to solve it. If you can have a look. Its easy but bit tricky. – Usman YousafZai Jan 05 '18 at 15:25
  • 1
    There are tidier ways, but to make it fit in this comment, `ii = r.argsort(); Pbest = Pbest.T[ii].T; r = r[ii]` – Eric Jan 05 '18 at 15:28
  • @Eric ... Its giving me the same error as it was giving before. setting an array element with a sequence. – Usman YousafZai Jan 05 '18 at 15:54
  • 1
    You should ask another question about simplifying `Ordenamiento_Burbuja`, showing what you tried, and what error you get. Make sure to remove any irrelevant code, like – Eric Jan 05 '18 at 15:55
  • @Eric ... will ask but I am not sure about that whether I made more mistakes or made it correct. I will edit the same question with ode libraries as the pso code above is working fine by without integrating the ode function. – Usman YousafZai Jan 05 '18 at 15:59
  • Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/162615/discussion-between-a4-page-and-eric). – Usman YousafZai Jan 05 '18 at 16:09

0 Answers0