-- encoding: utf-8 --
import torch_geometric
import argparse
from ast import arg
import torch
import math
import numpy as np
import pandas as pd
import torch.nn as nn
from torch.utils.data import Dataset, DataLoader
import copy
import os
import sys
import torch.utils.data as Data
#from optuna import TrialState
import optuna
import torch.nn.functional as F
#from get_data import get_mape
#from model import GAT_mlpModel
from itertools import chain
#from sklearn.preprocessing import MinMaxScaler
import torch
from scipy.interpolate import make_interp_spline
from torch import nn
import numpy as np
import matplotlib.pyplot as plt
from torch.utils.data import DataLoader
import optuna
from pandas import Series
from matplotlib import pyplot
from tqdm import tqdm
from torch.optim.lr_scheduler import StepLR
from tqdm import tqdm
#from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error,mean_absolute_percentage_error
import torch
from itertools import chain
from torch import nn
from torch_geometric.nn import GATConv, GCNConv, SAGEConv
from torch_geometric.data import Data #as graphData
from torch_geometric.loader import DataLoader
from torch_geometric.utils import to_undirected
#from get_data import nn_seq
#from args import args_parser
#from util import train, test, get_best_parameters
plt.rcParams[‘font.sans-serif’]=[‘SimHei’] #用来正常显示中文标签
plt.rcParams[‘axes.unicode_minus’]=False #用来正常显示负号 #有中文出现的情况,需要u’内容’
import scipy.sparse as sp
import os
os.environ[“KMP_DUPLICATE_LIB_OK”]=“TRUE”
path = os.path.dirname(os.path.realpath(file))
device = torch.device(“cuda” if torch.cuda.is_available() else “cpu”)
def gnn_args_parser():
parser = argparse.ArgumentParser()
parser.add_argument('--epochs', type=int, default=10, help='training epochs')
parser.add_argument('--input_size', type=int, default=27, help='input dimension')
parser.add_argument('--seq_len', type=int, default=16, help='seq len')
parser.add_argument('--output_size', type=int, default=16, help='output dimension')
parser.add_argument('--hidden_size', type=int, default=64, help='hidden size')
parser.add_argument('--lr', type=float, default=0.008, help='learning rate')#0.008
parser.add_argument('--batch_size', type=int, default=128, help='batch size')
parser.add_argument('--optimizer', type=str, default='adam', help='type of optimizer')
parser.add_argument('--device', default=torch.device("cuda" if torch.cuda.is_available() else "cpu"))
parser.add_argument('--weight_decay', type=float, default=1e-4, help='weight decay')
parser.add_argument('--step_size', type=int, default=150, help='step size')
parser.add_argument('--gamma', type=float, default=0.5, help='gamma')
args = parser.parse_args()
return args
def calc_corr(a, b):
s1 = Series(a)
s2 = Series(b)
return s1.corr(s2)
class MyDataset(Dataset):
def init(self, data):
self.data = data
def __getitem__(self, item):
return self.data[item]
def __len__(self):
return len(self.data)
def create_graph(num_nodes, data):
features = torch.randn((num_nodes, 256))
edge_index = [[], []]
# 计算相关系数
# data (x, num_nodes)
for i in range(num_nodes):
for j in range(i + 1, num_nodes):
x, y = data[:, i], data[:, j]
corr = calc_corr(x, y)
if corr >= 0.4:
edge_index[0].append(i)
edge_index[1].append(j)
edge_index = torch.LongTensor(edge_index)
graph = Data(x=features, edge_index=edge_index)
graph.edge_index = to_undirected(graph.edge_index, num_nodes=num_nodes)
# print(graph)
return graph
def adj2coo(adj):
# adj numpy
edge_index_temp = sp.coo_matrix(adj)
values = edge_index_temp.data
indices = np.vstack((edge_index_temp.row, edge_index_temp.col))
edge_index = torch.LongTensor(indices)
return edge_index
def nn_seq_gnn(num_nodes, seq_len, B, pred_step_size):
filename = ‘/data/data.csv’
filepath = os.path.dirname(os.path.realpath(file)) + filename
data = pd.read_csv(filepath, parse_dates=True, index_col=r’时间’)
data = data.fillna(‘0’)
data[“cos预测风向”]=np.cos(data[“预测风向(°)”].values.astype(‘float’))
data[“sin预测风向”]=np.sin(data[“预测风向(°)”].values.astype(‘float’))
data[“实际功率”] = data[“实际功率(MW)”]
#实际功率(MW),平均风速1(m/s),最大风速1(m/s),极大风速1(m/s),平均风向1(°),平均风速2(m/s),最大风速2(m/s),极大风速2(m/s),平均风向2(°),平均风速3(m/s),最大风速3(m/s),极大风速3(m/s),平均风向3(°),平均风速4(m/s),最大风速4(m/s),极大风速4(m/s),平均风向4(°),平均风速5(m/s),最大风速5(m/s),极大风速5(m/s),平均风向5(°),平均风速6(m/s),最大风速6(m/s),极大风速6(m/s),平均风向6(°),平均气温(℃),平均湿度(%),平均气压(hPa),平均空气密度(kg/m3),预测风速(m/s),预测风向(°),预测气温(℃),
data = data[[“实际功率”,“平均风速1(m/s)”,“最大风速1(m/s)”,“极大风速1(m/s)”,“平均风向1(°)”,“平均风速2(m/s)”,“最大风速2(m/s)”,“极大风速2(m/s)”,“平均风向2(°)”,“平均风速3(m/s)”,“最大风速3(m/s)”,“极大风速3(m/s)”,“平均风向3(°)”,“平均风速4(m/s)”,“最大风速4(m/s)”,“极大风速4(m/s)”,“平均风向4(°)”,“平均风速5(m/s)”,“最大风速5(m/s)”,“极大风速5(m/s)”,“平均风向5(°)”,“预测风速(m/s)”,“预测风向(°)”,“预测气温(℃)”]]
#data = data[[“实际功率”,“预测风速(m/s)”,“预测风向(°)”,“预测气温(℃)”,“cos预测风向”,“sin预测风向”]]#
train = data[:‘2024-01-01 00:15:00’]
val = data[‘2023-12-01 00:15:00’:‘2024-01-01 00:15:00’]
test = data[‘2023-12-31 16:30:00’:]
result=pd.DataFrame(columns=['实际功率', '真实值', '预测值'])
result['实际功率']=test['实际功率']['2024-01-01 00:15:00':]#[96+96-1:]
print(train.shape,val.shape,test.shape)
print(train.shape,test.shape)
num_nodes=24#3
#print('9'*100)
#print("train",train.shape)
graph = create_graph(num_nodes, train.values)
def process(dataset, batch_size, step_size, shuffle):
dataset = dataset.values.tolist()
graphs = []
ys=[]
for i in tqdm(range(0, len(dataset) - seq_len - pred_step_size+1, step_size)):
train_seq = []
for j in range(i, i + seq_len):
x = []
for c in range(len(dataset[0])): # 前24个时刻的所有变量
x.append(dataset[j][c])
train_seq.append(x)
# 下几个时刻的所有变量
train_labels = []
for j in range(len(dataset[0])):
train_label = []
for k in range(i + seq_len, i + seq_len + pred_step_size):
train_label.append(dataset[k][j])
train_labels.append(train_label)
# tensor
#print('o'*100)
#print('train_seq',train_seq)
train_seq=np.array(train_seq).astype("float")
train_labels=np.array(train_labels).astype("float")
train_seq = torch.FloatTensor(train_seq)
train_labels = torch.FloatTensor(train_labels)
# print(train_seq.shape, train_labels.shape) # 24 13, 13 1
# 此处可利用train_seq创建动态的邻接矩阵
temp = Data(x=train_seq.T, edge_index=graph.edge_index, y=train_labels)
# print(temp)
y_temp=Data(y=train_labels)
graphs.append(temp)
ys.append(y_temp)
loader = torch_geometric.loader.DataLoader(graphs, batch_size=batch_size,
shuffle=shuffle, drop_last=False)
return loader
Dtr = process(train, B, step_size=1, shuffle=False)
Val = process(val, B, step_size=1, shuffle=False)
Dte = process(test, B, step_size=1, shuffle=False)
return graph, Dtr, Val, Dte,result
class GAT(torch.nn.Module):
def init(self, in_feats, h_feats, out_feats):
super(GAT, self).init()
self.conv1 = GATConv(in_feats, h_feats, heads=4, concat=False)
self.conv2 = GATConv(h_feats, out_feats, heads=4, concat=False)
def forward(self, x, edge_index, edge_weight=None):
# 24 128 / 2 118
x = F.elu(self.conv1(x, edge_index))
x = self.conv2(x, edge_index)
return x
class GNN_MLP(nn.Module):
def init(self, args):
super(GNN_MLP, self).init()
self.args = args
self.out_feats = 128
self.gat = GAT(in_feats=args.seq_len, h_feats=100, out_feats=self.out_feats)
self.fc = nn.Sequential(
nn.Linear(128, 64),
nn.ReLU(),
nn.Linear(64, args.output_size)
)
self.fcs = nn.ModuleList()
for k in range(args.input_size):
self.fcs.append(nn.Sequential(
nn.Linear(128, 64),
nn.ReLU(),
nn.Linear(64, args.output_size)
))
def forward(self, data):
# Data(x=[13, 24], edge_index=[2, 32], y=[13, 1])
# DataBatch(x=[6656, 24], edge_index=[2, 16384], y=[6656, 1], batch=[6656], ptr=[513])
# output(13, 512, 1) y(512, 13, 1)
x, edge_index, batch = data.x, data.edge_index, data.batch
batch_size = torch.max(batch).item() + 1
x = self.gat(x, edge_index) # 6656 128 = 512 * (13, 128) # y = 6656 1 = 512 * (13 1)
batch_list = batch.cpu().numpy()
# print(batch_list)
# split
ys = [[] for k in range(batch_size)]
xs = [[] for k in range(batch_size)]
for k in range(x.shape[0]):
xs[batch_list[k]].append(x[k, :])
ys[batch_list[k]].append(data.y[k, :])
xs = [torch.stack(x, dim=0) for x in xs]
x = torch.stack(xs, dim=0)
#print('p'*100)
#print(x.shape) # 512 13 128 / 512 13 1
# output(13, 512, 1) y(512, 13, 1)
preds = []
x = x.permute(1, 0, 2) # 13 512 128
for idx, fc in enumerate(self.fcs):
preds.append(fc(x[0, :, :]))
break
pred = torch.stack(preds, dim=0)
ys = [torch.stack(x, dim=0) for x in ys]
y = torch.stack(ys, dim=0)
'''
ys = [[] for k in range(batch_size)]
for k in range(x.shape[0]):
xs[batch_list[k]].append(x[k, :])
ys[batch_list[k]].append(data.y[k, :])
xs = [torch.stack(x, dim=0) for x in xs]
ys = [torch.stack(x, dim=0) for x in ys]
x = torch.stack(xs, dim=0)
y = torch.stack(ys, dim=0)
'''
y = y[:, 0, :]
y = y.unsqueeze(1)
#print(pred.shape,y.shape)
return pred,y
def print_history(train_loss_,val_loss_,name):
plt.plot(train_loss_,color=‘Orange’,label=“train loss”)
#print(“history.history[‘loss’]”,len(history.history[‘loss’]))
plt.plot(val_loss_,color=‘b’,label=“validation loss”)
plt.title(‘train_validation loss’)
plt.xlabel(‘epoch’)
plt.ylabel(‘loss’)
plt.legend()
#plt.show()
plt.savefig(path+r’/loss/‘+name+’.png’)
plt.close()
def get_val_loss2(args, model, Val):
model.eval()
loss_function = nn.MSELoss().to(args.device)
val_loss = []
for graph in Val:
graph = graph.to(device)
preds, labels = model(graph)
total_loss = 0
#for k in range(args.input_size):
total_loss = total_loss + loss_function(preds[0, :, :], labels[:,0, :])
total_loss /= preds.shape[0]
val_loss.append(total_loss.item())
return np.mean(val_loss)
def train2(args, Dtr, Val):
model_type = ‘gnn2’
model = GNN_MLP(args).to(device)
name=str(‘GNN2’)
loss_function = nn.MSELoss().to(device)
if args.optimizer == 'adam':
optimizer = torch.optim.Adam(model.parameters(), lr=args.lr,
weight_decay=args.weight_decay)
else:
optimizer = torch.optim.SGD(model.parameters(), lr=args.lr,
momentum=0.9, weight_decay=args.weight_decay)
scheduler = StepLR(optimizer, step_size=args.step_size, gamma=args.gamma)
# training
min_epochs = 2
best_model = model
min_val_loss = 5
train_loss_,val_loss_=[],[]
for epoch in tqdm(range(args.epochs)):
train_loss = []
for graph in Dtr:
graph = graph.to(device)
preds,labels = model(graph)
total_loss = 0
#for k in range(args.input_size):
total_loss = total_loss + loss_function(preds[0, :, :], labels[:, 0, :])
total_loss = total_loss / preds.shape[0]
optimizer.zero_grad()
total_loss.backward()
optimizer.step()
train_loss.append(total_loss.item())
scheduler.step()
# validation
val_loss = get_val_loss2(args, model, Val)
if epoch + 1 >= min_epochs and val_loss < min_val_loss:
min_val_loss = val_loss
best_model = copy.deepcopy(model)
state = {'model': best_model.state_dict()}
torch.save(state, path + '/models/' + model_type + '.pkl')
print('epoch {:03d} train_loss {:.8f} val_loss {:.8f}'.format(epoch, np.mean(train_loss), val_loss))
train_loss_.append(np.mean(train_loss))
val_loss_.append(val_loss)
model.train()
print_history(train_loss_,val_loss_,name)
#state = {'model': best_model.state_dict()}
torch.save(best_model, path + '/models/' + model_type + '.pkl')
@torch.no_grad()
def test2(args, Dte, Cap,result):
print(‘loading models…’)
model_type = ‘gnn2’
model = GNN_MLP(args).to(device)
model=torch.load(path + ‘/models/’ + model_type + ‘.pkl’).to(device)#, map_location=torch.device(‘cpu’))
model.eval()
print(‘predicting…’)
targets_list=[]
preds_list=[]
for graph in tqdm(Dte):
graph = graph.to(device)
preds, targets = model(graph)
targets = torch.from_numpy(np.array(targets.data.tolist())) # (batch_size, n_outputs, pred_step_size)
targets=targets.squeeze(dim=1)
preds=preds.squeeze(dim=0)
targets_list.append(targets)
preds_list.append(preds)
#print(targets.shape,preds.shape)
‘’’
for i in range(args.input_size):
target = targets[:, i, :]
target = list(chain.from_iterable(target))
ys[i].extend(target)
break
for i in range(_pred.shape[0]):
pred = _pred[i]
pred = list(chain.from_iterable(pred.data.tolist()))
preds[i].extend(pred)
break
‘’’
# ys, preds = [np.array(y) for y in ys], [np.array(pred) for pred in preds]
#targets=targets.squeeze(dim=1)
#preds=preds.squeeze(dim=0)
targets=torch.cat(targets_list, dim=0)
preds=torch.cat(preds_list,dim=0)
#print(targets.shape,preds.shape)
#ys = scaler.inverse_transform(ys).T
#preds = scaler.inverse_transform(preds).T
preds = preds.detach().cpu().numpy()
ys = targets.numpy()
#print(ys.shape,preds.shape)
#y = scaler1.inverse_transform(ys)#.T
#pred = scaler1.inverse_transform(preds)#.T
y,pred=ys[:,-1],preds[:,-1]
for j in range(len(pred)):
pred[j] = np.round(pred[j], 3)
if pred[j] < 0:
pred[j] = float(0)
if pred[j]>Cap:
pred[j]=Cap
print(y.shape,pred.shape)
#print(type(y),type(pred))
print(‘mse:’, get_mse(y, pred))
print(‘rmse:’, get_rmse(y, pred))
print(‘mae:’, get_mae(y, pred))
print(‘mape:’, get_mape(y, pred))
print(‘r2score’,get_r2(y,pred))
result[‘真实值’]=y#[2192-2014:]
result[‘预测值’]=pred#[2192-2014:]
plot_test(y, pred)
result.to_csv(path + ‘/result/’ + model_type + ‘-result.csv’, sep=‘,’)
mse=get_mse(y, pred)
rmse=get_rmse(y, pred)
mae=get_mae(y, pred)
mape= get_mape(y, pred)
r2score=get_r2(y,pred)
df = {“mse”: mse, “rmse”: rmse,“mae”: mae, “mape”: mape,“r2score”:r2score}
df = pd.DataFrame(list(df.items()))
df.to_csv(path + ‘/result_score/’ + model_type + ‘-result_score.csv’)
@torch.no_grad()
def test16(args, Dte, Cap,result):#,scaler1):# scaler,
print(‘loading models…’)
model_type = ‘gnn2’
model = GNN_MLP(args).to(device)
model=torch.load(path + ‘/models/’ + model_type + ‘.pkl’).to(device)#, map_location=torch.device(‘cpu’))
model.eval()
print(‘predicting…’)
targets_list=[]
preds_list=[]
for graph in tqdm(Dte):
graph = graph.to(device)
preds, targets = model(graph)
targets = torch.from_numpy(np.array(targets.data.tolist())) # (batch_size, n_outputs, pred_step_size)
targets=targets.squeeze(dim=1)
preds=preds.squeeze(dim=0)
targets_list.append(targets)
preds_list.append(preds)
#print(targets.shape,preds.shape)
‘’’
for i in range(args.input_size):
target = targets[:, i, :]
target = list(chain.from_iterable(target))
ys[i].extend(target)
break
for i in range(_pred.shape[0]):
pred = _pred[i]
pred = list(chain.from_iterable(pred.data.tolist()))
preds[i].extend(pred)
break
‘’’
# ys, preds = [np.array(y) for y in ys], [np.array(pred) for pred in preds]
#targets=targets.squeeze(dim=1)
#preds=preds.squeeze(dim=0)
targets=torch.cat(targets_list, dim=0)
preds=torch.cat(preds_list,dim=0)
print(targets.shape,preds.shape)
#ys = scaler.inverse_transform(ys).T
#preds = scaler.inverse_transform(preds).T
preds = preds.detach().cpu().numpy()
ys = targets.numpy()
#print(ys.shape,preds.shape)
#y = scaler1.inverse_transform(ys)#.T
#pred = scaler1.inverse_transform(preds)#.T
y,pred=ys[:,-1],preds[:,-1]
print(y.shape,pred.shape)
#print(type(y),type(pred))
print(‘mse:’, get_mse(y, pred))
print(‘rmse:’, get_rmse(y, pred))
print(‘mae:’, get_mae(y, pred))
print(‘mape:’, get_mape(y, pred))
print(‘r2score’,get_r2(y,pred))
result[‘真实值’]=y#[2192-2014:]
result[‘预测值’]=pred#[2192-2014:]
plot_test(y, pred)
def correction(jj):
for j in range(len(preds[:,jj])):
preds[:,jj][j] = np.round(preds[:,jj][j], 3)
if preds[:,jj][j] < 0:
preds[:,jj][j] = float(0)
if preds[:,jj][j]>Cap:
preds[:,jj][j]=Cap
for j in range(16):
correction(j)
result['真实值']=y
for i in range(16):
result['预测值'+str(i)]=preds[:,i]
result.to_csv(path + '/result/' + model_type + '16个点-result.csv', sep=',')
mse=get_mse(y, pred)
rmse=get_rmse(y, pred)
mae=get_mae(y, pred)
mape= get_mape(y, pred)
r2score=get_r2(y,pred)
df = {"mse": mse, "rmse": rmse,"mae": mae, "mape": mape,"r2score":r2score}
df = pd.DataFrame(list(df.items()))
df.to_csv(path + '/result_score/' + model_type + '取16个点最后一个-result_score.csv')
def plot_test(y, pred):
# plot
plt.plot(y, color=‘blue’, label=‘true value’)
plt.plot(pred, color=‘red’, label=‘pred value’)
plt.title(‘GNN的预测结果’)
#plt.grid(True)
plt.legend(loc=‘upper center’, ncol=6)
#plt.show()
plt.savefig(path+r’/pictures/gnn2.png’)
def get_mape(y_true, y_pred):
“”"
计算平均绝对百分比误差 (Mean Absolute Percentage Error)
“”"
# 防止除以零和负数开根号的情况
y_true = np.where(y_true == 0, 1e-6, y_true) # 对于0值,用一个极小正数替代以避免除法错误
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
def get_r2(y_true, y_pred):
“”"
计算决定系数 R²
“”"
rss = np.sum((y_true - y_pred) ** 2) # 剩余平方和
tss = np.sum((y_true - np.mean(y_true)) ** 2) # 总平方和
return 1 - (rss / tss)
def get_mae(y_true, y_pred):
“”"
计算平均绝对误差 (Mean Absolute Error)
“”"
return np.mean(np.abs(y_true - y_pred))
def get_mse(y_true, y_pred):
“”"
计算均方误差 (Mean Squared Error)
“”"
return np.mean((y_true - y_pred) ** 2)
def get_rmse(y, pred):
return np.sqrt(get_mse(y, pred))
def get_best_parameters(args, Dtr, Val):
def objective(trial):
model = GNN_MLP(args).to(args.device)
loss_function = nn.MSELoss().to(args.device)
optimizer = trial.suggest_categorical(‘optimizer’,
[torch.optim.SGD,
torch.optim.RMSprop,
torch.optim.Adam])(
model.parameters(), lr=trial.suggest_loguniform(‘lr’, 5e-4, 1e-2))
print(‘training…’)
epochs = 10
val_loss = 0
for epoch in range(epochs):
train_loss = []
for batch_idx, graph in enumerate(Dtr, 0):
graph = graph.to(args.device)
optimizer.zero_grad()
y_pred = model(graph).to(args.device)
target=graph.y.to(args.device)
#print(y_pred.shape,target.shape)
#target=target.reshape(-1,1)
loss = loss_function(y_pred, target)
train_loss.append(loss.item())
loss.backward()
optimizer.step()
# validation
val_loss = get_val_loss2(args, model, Val)
print('epoch {:03d} train_loss {:.8f} val_loss {:.8f}'.format(epoch, np.mean(train_loss), val_loss))
model.train()
return val_loss
sampler = optuna.samplers.TPESampler()
study = optuna.create_study(sampler=sampler, direction='minimize')
study.optimize(func=objective, n_trials=5)
pruned_trials = study.get_trials(deepcopy=False,
states=tuple([optuna.trial.TrialState.PRUNED]))
complete_trials = study.get_trials(deepcopy=False,
states=tuple([optuna.trial.TrialState.COMPLETE]))
best_trial = study.best_trial
print('val_loss = ', best_trial.value)
for key, value in best_trial.params.items():
print("{}: {}".format(key, value))
def main():
Cap=300
args = gnn_args_parser()
graph, Dtr, Val, Dte,result= nn_seq_gnn(args.input_size,args.seq_len,args.batch_size,args.output_size)
#get_best_parameters(args, Dtr, Val)
import time
T1 = time.time()
train2(args, Dtr, Val)#, PATH
print(‘train用时’,time.time()-T1)
T2=time.time()
test2(args, Dte, Cap,result)
test16(args, Dte, Cap,result)
print(‘test用时’,time.time()-T2)
if name == ‘main’:
main()