粒子群優化BP神經網路初始權值（python實作）

粒子群優化BP神經網路初始權值（python實作）

網上看了一些資料，但都是用matlab寫的，（還要用csdn會員積分下載）自己不太會用matlab，就試著用python寫了段小程式實作，資料用的是sklearn中的波士頓房價資料集，神經網路部分是用tensorflow2，（本來想用sklearn，但不知道怎么把網路權值提取出來）
思路挺簡單的，就是把BP網路誤差作為粒子群優化的目標，每個粒子對應網路初始權重，一輪一輪的迭代，（先是在一篇建模論文里看到了這個想法，雖然好像沒啥實用價值）
大一，只會捏泥巴，各位看官就全當茶余飯后圖個樂吧，

#import tensorflow.compat.v1 as tf
#tf.compat.v1.disable_v2_behavior()
#import tensorflow as tf
#第一次用上面的陳述句跑的時候還好好的，再跑就報錯了
import tensorflow as tf
tf = tf.compat.v1
tf.disable_v2_behavior()


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler


def function(x1,y1,x2,y2,W):#  W是神經網路的W權重，根據這個權重設定神經網路




    #定義激活函式
    activation_function=tf.nn.relu
    #輸入輸出資料集
    xs=tf.placeholder(tf.float32,[None,None])
    ys=tf.placeholder(tf.float32,[None,None])




    #設計bp神經網路，三層，13,3,1
    weights_1=tf.Variable(W[0,:,:],tf.float32)
    biases_1=tf.Variable(tf.zeros([1,3])+0.1,tf.float32)
    wx_plus_b_1=tf.matmul( xs, tf.cast(weights_1,tf.float32))+biases_1
    outputs_1=activation_function(wx_plus_b_1)


    weights_2=tf.Variable(W[1,0:3,:],tf.float32)
    biases_2=tf.Variable(tf.zeros([1,3])+0.1,tf.float32)
    wx_plus_b_2=tf.matmul(outputs_1 , tf.cast(weights_2,tf.float32))+biases_2
    outputs_2=activation_function(wx_plus_b_2)


    w3=W[2,0:3,0].reshape(3,1)
    weights_3=tf.Variable(w3,tf.float32)
    biases_3=tf.Variable(0.1,tf.float32)

    wx_plus_b_3=tf.matmul(outputs_2,tf.cast(weights_3,tf.float32))+biases_3



    #預測輸出結果
    prediction=wx_plus_b_3      #看來這里的資料就用行向量來輸入輸出

    #定義損失函式
    loss=tf.reduce_mean(tf.reduce_sum(tf.square(y1-prediction),reduction_indices=[1]))

    #梯度下降法訓練
    train_step=tf.train.GradientDescentOptimizer(0.1).minimize(loss)

    #初始化變數
    init=tf.global_variables_initializer()

    #執行會話，開始訓練模型
    print("開始")
    with tf.Session() as sess:
        sess.run(init)
        for i in range (1000):
            sess.run(train_step,feed_dict={ xs:x1  , ys:y1 })

#為什么損失函式喂入x2,y2就不行？QAQ
        end_loss=sess.run(loss,feed_dict={xs:x1,ys:y1})
        print(end_loss)
# print(sess.run(prediction,feed_dict={xs:x2}))
        print("結束")
    return end_loss


#匯入資料集
data=load_boston()
data_pd=pd.DataFrame(data.data,columns=data.feature_names)
data_pd["price"]=data.target




#dataframe匯入numpy
x=np.array(data_pd.loc[:,'CRIM':'LSTAT'])

y=np.array(data_pd.loc[:,'price'])

y.shape=(506,1)
#訓練集測驗集
x_train,x_test,y_train,y_test=train_test_split(x,y , test_size=0.1 )
#資料標準化
SC=StandardScaler()
x_train=SC.fit_transform(x_train)
y_train=SC.fit_transform(y_train)
x_test=SC.fit_transform(x_test)
y_test=SC.fit_transform(y_test)



#粒子數量num
num = 3

#粒子位置矩陣的形狀
num_x = 3
num_y = 13
num_z = 3

#p為粒子位置矩陣，初始化為標準正態分布
p = np.random.randn(num,num_x,num_y,num_z)

#初始化粒子速度,以標準正態分布隨機初始化
v = np.random.randn(num,num_x,num_y,num_z)

#個體最佳位置
good_p = np.array(p, copy=True)

#全域最佳位置
best_p = np.zeros((num_x, num_y, num_z))

#每次粒子移動后所計算出新的目標函式值
new_y = np.zeros(num)

#粒子個體歷史最優值
good_y = np.zeros(num)

#粒子群體歷史最優值
best_y = 0

#計算出初始粒子群的目標函式值
for i in range(num):
    good_y[i] = function(x_train, y_train, x_test, y_test, p[i, :, :, :])

#目標函式回傳值是誤差，那么最小的就是最優的
best_y = min(good_y)

#確定初始時最優位置
best_p = p[np.argmin(good_y), :, :, :]

#設定最大迭代次數
max_iter = 10

#開始迭代
for i in range(max_iter):

    #速度更新公式
    v = random.random() * v + 2.4 * random.random() * (best_p - p) + 1.7 * random.random() * ( good_p - p )

    #粒子位置更新
    p = p + v

    #計算每個粒子到達新位置后所得到的目標函式值
    for i in range(num):
        new_y[i] = function(x_train, y_train, x_test, y_test, p[i, :, :, :])

    #更新全域最優
    if min(new_y) < best_y:
        best_y = min(new_y)
        best_p = p[np.argmin(new_y), :, :, :]

    #更新個體歷史最優
    for i in range(num):
        if new_y[i] < good_y[i]:
            good_y[i] = new_y[i]
            good_p[i, :, :, :] = p[i, :, :, :]  # 當對切片修改時，原始numpy資料也修改


print("結束")
print('目標函式最優值：',best_y)
print('此時的粒子位置：',best_p)