本人很小白，SA和BPNN的原理都比较清楚，但是搜了一下发现没有用SA优化BPNN的Python代码，毕业论文急用
附上我的BPNN和SA分别的代码
BPNN：
# -*- coding: utf-8 -*-
from __future__ import print_function
import theano
import theano.tensor as T
from theano.ifelse import ifelse
import numpy as np
import pandas as pd
import os
from sklearn.metrics import mean_squared_error, r2_score,mean_absolute_error
import warnings
warnings.filterwarnings('ignore')
import time
start=time.clock()
#定义神经网络模型：反向传播神经网络；输入层/隐藏层/输出层；
class Layer(object):
def __init__(self, inputs, in_size, out_size, activation_function=None):
self.W = theano.shared(np.random.normal(0, 1, (in_size, out_size)))
self.b = theano.shared(np.zeros((out_size, )) + 0.1)
self.Wx_plus_b = T.dot(inputs, self.W) + self.b
self.activation_function = activation_function
if activation_function is None:
self.outputs = self.Wx_plus_b
else:
self.outputs = self.activation_function(self.Wx_plus_b)
# import data
os.chdir(r'C:\Users\Songuoer\Desktop\不知道\订单数据')
df = pd.read_excel('净味防霉随机测试数据2018.xlsx',sheet_name='data')
df['lnq'] = np.log(df['q'])
df['q_1'] = df['q'].shift(1)
df['lnq_1'] = df['lnq'].shift(1)
df.dropna(inplace=True)
# 数据标准化
def max_min(data):
data0 = [(x - min(data))/(max(data) - min(data)) for x in data]
return data0
# 标准化后的data
data = df[['ym']]
data['p'] = max_min(df['p'])
data['discount'] = max_min(df['discount'])
data['season'] = max_min(df['season'])
data['q'] = max_min(df['q'])
data['q_1'] = max_min(df['q_1'])
#定义输入和输出矩阵：
feature_cols = ['q_1','p','discount','season']
# 测试集
x_test = np.array(data[feature_cols][:40])
y_test = np.array(data['q'])[:40][:, np.newaxis]
# 训练集
x_pre = np.array(data[feature_cols][40:])
y_pre = np.array(data['q'])[40:][:, np.newaxis]
# 实际销量
q_test = np.array(df['q'][40:],dtype=float)
#定义神经网络内的bias调整偏差值
noise = np.random.normal(0, 0.35, x_test.shape)
# determine the inputs dtype
x = T.dmatrix("x")
y = T.dmatrix("y")
#定义神经网络的隐藏层和输出层：
# add layers
hidden = 8
l1 = Layer(x, 4, hidden, T.nnet.sigmoid)
l2 = Layer(l1.outputs,hidden, 1, None)
# compute the cost
cost = T.mean(T.square(l2.outputs - y))
'''
#这里说明一下希望把cost函数改成这种形式
def cost(l2outputs,y):
y_gap = l2outputs-y
yh = y_gap[y_gap>0]*2
yw = -y_gap[y_gap<0]*10
tc = T.sum(yh)+T.sum(yw)
return tc
'''
# compute the gradients
gW1, gb1, gW2, gb2 = T.grad(cost, [l1.W, l1.b, l2.W, l2.b])
# apply gradient descent
learning_rate = 0.001
train = theano.function(
inputs=[x, y],
outputs=cost,
updates=[(l1.W, l1.W - learning_rate * gW1),
(l1.b, l1.b - learning_rate * gb1),
(l2.W, l2.W - learning_rate * gW2),
(l2.b, l2.b - learning_rate * gb2)])
# prediction
predict = theano.function(inputs=[x], outputs=l2.outputs)
# BPNN的参数
paramters = theano.function(inputs=[x], outputs=[l1.W, l1.b, l2.W, l2.b]
,on_unused_input='ignore')
# 迭代预测
total_times = 30000
for i in range(total_times):
# training
err = train(x_test, y_test)
if i % 6000 == 0:
print('%i times error:'%i,err)
prediction_value = predict(x_pre)
par = paramters(x_pre)
# 将预测值还原
def nor_q(q):
minq,maxq = min(df['q']),max(df['q'])
q_o = [((maxq-minq)*x+minq) for x in q]
return q_o
q_pre = np.array(nor_q(prediction_value)).reshape(7,)
final_par_z = par
# 定义MAPE
def mape(y_true, y_pred):
return np.mean(np.abs((y_pred - y_true) / y_true)) * 100
print('Variance score: %.3f' % r2_score(q_test, q_pre))
print("RMSE: %.2f" % np.sqrt(mean_squared_error(q_test, q_pre)))
print('MAE: %.3f'%mean_absolute_error(q_test,q_pre))
print('MAPE: %.3f'%mape(q_test, q_pre))
result = pd.DataFrame({'实际销量':q_test,'预测销量':q_pre},
index=pd.period_range('6/1/2018','12/1/2018',freq = 'M'))
end=time.clock()
print('Running time: %s Seconds'%(end-start))
print('finished!')
SA：
import numpy as np
import pandas as pd
import math
import os
import warnings
warnings.filterwarnings('ignore')
# 不发出警告
import time
start=time.clock()
os.chdir(r'C:\Users\Songuoer\Desktop\不知道\订单数据')
df = pd.read_excel('真时丽随机测试数据2018.xlsx',sheet_name='data')
data = df[40:]
data['lnq'] = np.log(data['q'])
data['lnq_1'] = data['lnq'].shift()
h = 2
w = 12
beta1 = [-12.5639,0.2933,0.184,-0.0476,0.8487] #真时丽
beta2 = [7.8332,0.2464,0.0584,-0.1038,0.0178] #家净丽
beta3 = [20.604908,0.010558,-0.165661,0.016112] #净味防霉，没有lnq_1
beta = beta1
#define aim function
def aimFunction(beta):
data['lnq_pre'] = beta[0]+beta[1]*data['lnq_1']+beta[2]*data['p']+beta[3]*data['discount']+beta[4]*data['season']
# data['lnq_pre'] = beta[0]+beta[1]*data['p']+beta[2]*data['discount']+beta[3]*data['season']
data['lnq_pre'][40] = data['lnq'][40]
data['OUL'] = np.exp(data['lnq_pre']).astype(np.int)
data['inventory'] = data['OUL']-data['q']
data['short'] = data['q']-data['OUL']
data['inventory'][data['inventory']<0] = 0
data['short'][data['short']<0] = 0
data['C_inventory'] = h*data['inventory']
data['C_short'] = w*data['short']
data['Cost'] = data['C_inventory']+data['C_short']
TC = data['Cost'].iloc[1:].sum()
return TC
TC0 = aimFunction(beta)#initiate result
print('initiate result:%i'%TC0)
#模拟退火算法
T=1500 #initiate temperature
Tmin=10 #minimum value of terperature
betaNew = [0,0,0,0,0]
k=200 #times of internal circulation
t=0 #time
a = 0.2
while T>=Tmin:
for i in range(k):
#calculate y
TC=aimFunction(beta)
#generate a new x in the neighboorhood of x by transform function
for i in range(5):
if beta[i]>0:
betaNew[i] = np.random.uniform(low=(1-a)*beta[i],high=(1+a)*beta[i])
else:
betaNew[i] = np.random.uniform(low=(1+a)*beta[i],high=(1-a)*beta[i])
TCNew=aimFunction(betaNew)
if TCNew-TC<0:
beta=betaNew.copy()
else:
#metropolis principle
p=math.exp(-(TCNew-TC)/T)
r=np.random.uniform(low=0,high=1)
if r<p:
beta=betaNew.copy()
t+=1
if t%20 == 0:
print ('Intervel: %i,Total Cost:%i'%(t,aimFunction(beta)))
T=1000/(1+t)
print (beta,aimFunction(beta))
end=time.clock()
print('Running time: %s Seconds'%(end-start))
希望能用模拟退火算法来优化BPNN的参数从而求最小总成本
谢谢各位