This commit is contained in:
chenxiaodong 2024-07-30 09:05:32 +08:00
parent a6232f24c8
commit 23c1ad592d
19 changed files with 3183 additions and 6737 deletions

147
7.24.py Normal file
View File

@ -0,0 +1,147 @@
import numpy as np
import pandas as pd
import time
import os
import json
from concurrent.futures import ThreadPoolExecutor
from deap import base, creator, tools, algorithms
def fitness_numpy(individual, price, load, temperature, irradiance, wind_speed, prev_soc, prev_Pg1, prev_Pg2, prev_Pg3):
individual = np.array(individual)
price = np.array(price)
load = np.array(load)
temperature = np.array(temperature)
irradiance = np.array(irradiance)
wind_speed = np.array(wind_speed)
Ac, Ag1, Ag2, Ag3, Av = individual[:5]
soc = np.clip(prev_soc + 0.2 * Ac * 0.9, 0.2, 0.8)
Pg1 = np.clip(prev_Pg1 + 100 * Ag1, 0, 150)
Pg2 = np.clip(prev_Pg2 + 100 * Ag2, 0, 375)
Pg3 = np.clip(prev_Pg3 + 200 * Ag3, 0, 500)
Pso = np.clip((0.2 * irradiance + 0.05 * temperature - 9.25) * (1 + Av), 0, None)
Pw = np.where((wind_speed >= 3) & (wind_speed < 8), wind_speed ** 3 * 172.2625 / 1000,
np.where((wind_speed >= 8) & (wind_speed < 12), 64 * 172.2625 / 125, 0))
P = Ac + Pg1 + Pg2 + Pg3 + Pso + Pw
Ee = np.where(P >= load, P - load, 0)
Es = np.where(P < load, load - P, 0)
Cb = 0.01 * Ac + 0.1 * soc
Cg1 = 0.0034 * Pg1 ** 2 + 3 * Pg1 + 30
Cg2 = 0.001 * Pg2 ** 2 + 10 * Pg2 + 40
Cg3 = 0.001 * Pg3 ** 2 + 15 * Pg3 + 70
Cs = 0.01 * Pso
Cw = 0.01 * Pw
Rs = 0.5 * price * Ee
Cp = price * Es
Pe = np.where(Ee > 100, (Ee - 100) * 50, 0)
Ps = np.where(Es > 100, (Es - 100) * 50, 0)
total_cost = np.sum(Cb + Cg1 + Cg2 + Cg3 + Cs + Cw - Rs + Cp + Pe + Ps)
reward = -total_cost / 1000
return (reward,)
def check_bounds(func):
def wrapper(*args, **kwargs):
offspring = func(*args, **kwargs)
for individual in offspring:
for i in range(len(individual)):
individual[i] = np.clip(individual[i], -1, 1)
return offspring
return wrapper
def main():
data = pd.read_csv('./data.csv')
price = data['price'].values
load = data['load'].values
temperature = data['temperature'].values
irradiance = data['irradiance'].values
wind_speed = data['wind_speed'].values
creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("Individual", list, fitness=creator.FitnessMax)
toolbox = base.Toolbox()
toolbox.register("attr_float", np.random.uniform, -1, 1)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, n=5)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("mate", check_bounds(tools.cxBlend), alpha=0.5)
toolbox.register("mutate", check_bounds(tools.mutGaussian), mu=0, sigma=0.1, indpb=0.2)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("evaluate", fitness_numpy)
population = toolbox.population(n=500)
prev_soc, prev_Pg1, prev_Pg2, prev_Pg3 = 0.4, 0.0, 0.0, 0.0
decision_values = []
for hour in range(8760):
start = time.time()
current_price = price[hour]
current_load = load[hour]
current_temperature = temperature[hour]
current_irradiance = irradiance[hour]
current_wind_speed = wind_speed[hour]
for gen in range(500):
offspring = toolbox.select(population, len(population))
offspring = list(map(toolbox.clone, offspring))
for child1, child2 in zip(offspring[::2], offspring[1::2]):
if np.random.rand() < 0.7:
toolbox.mate(child1, child2)
del child1.fitness.values
del child2.fitness.values
for mutant in offspring:
if np.random.rand() < 0.2:
toolbox.mutate(mutant)
del mutant.fitness.values
invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
with ThreadPoolExecutor() as executor:
futures = []
batch_size = 500
for i in range(0, len(invalid_ind), batch_size):
batch = invalid_ind[i:i + batch_size]
for ind in batch:
futures.append(
executor.submit(toolbox.evaluate, ind, current_price, current_load, current_temperature,
current_irradiance, current_wind_speed, prev_soc, prev_Pg1, prev_Pg2,
prev_Pg3))
for future in futures:
fitness = future.result()
for ind in invalid_ind:
if not ind.fitness.valid:
ind.fitness.values = fitness
population[:] = offspring
end = time.time()
best_ind = tools.selBest(population, 1)[0]
print(f'Best individual at hour {hour + 1}: {best_ind}')
print(f'Fitness: {best_ind.fitness.values}, using {end - start}s')
decision_values.append({
'Ac': best_ind[0],
'Ag1': best_ind[1],
'Ag2': best_ind[2],
'Ag3': best_ind[3],
'Av': best_ind[4]
})
prev_soc, prev_Pg1, prev_Pg2, prev_Pg3 = best_ind[0], best_ind[1], best_ind[2], best_ind[3]
with open('decision_values.json', 'w') as f:
json.dump(decision_values, f)
if __name__ == "__main__":
main()

File diff suppressed because it is too large Load Diff

Before

Width:  |  Height:  |  Size: 145 KiB

File diff suppressed because it is too large Load Diff

Before

Width:  |  Height:  |  Size: 126 KiB

After

Width:  |  Height:  |  Size: 127 KiB

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

2
PPO.py
View File

@ -385,7 +385,7 @@ if __name__ == '__main__':
from plotDRL import PlotArgs, make_dir, plot_evaluation_information, plot_optimization_result
plot_args = PlotArgs()
plot_args.feature_change = 'gae'
plot_args.feature_change = 'gae_solar'
args.cwd = agent_name
plot_dir = make_dir(args.cwd, plot_args.feature_change)
plot_optimization_result(base_result, plot_dir)

56
Qwen_max.py Normal file
View File

@ -0,0 +1,56 @@
from openai import OpenAI
import os
def get_response():
client = OpenAI(
api_key="sk-f01744b2801344b1a72f89ec7e290cad",
base_url="https://dashscope.aliyuncs.com/compatible-mode/v1", # 填写DashScope服务的base_url
)
completion = client.chat.completions.create(
model="qwen-turbo",
messages=[
{'role': 'system',
'content':
'你是一名擅长调试和控制智慧住宅能源系统的专家,需要根据当前环境状态给出最优的动作决策。'
'能源系统结构包括太阳能发电站、风能发电站、发电机、电池模块和住宅电网连接。'},
{'role': 'user',
'content':
'一、当前状态包含9个值'
'Sc(电池荷电状态)初始为0.4范围0.2到0.8Sg1(发电机1功率状态)、Sg2(发电机2功率状态)和Sg3(发电机3功率状态)初始功率都为0'
'Sp(电价状态)、Sl(住宅负载状态)、St(温度状态)、Si(太阳辐射度状态)和Sw(风速状态)都由数据提供。'
'其中Sc、Sg1、Sg2和Sg3的值都是在上一时刻基础上进行变化的始终不小于0Sg1、Sg2和Sg3为0代表各发动机关闭。'
'二、当前动作包含5个值'
'x表示当前时间的取值, 其范围从-1到1电池总容量为500初始电池容量为500*0.4=200Ag1、Ag2和Ag3分别表示相对于各自最大瞬时功率的变化值。'
'Ac(电池充放电动作)电池容量变化100*xAg1(发电机1功率动作)功率变化100*xAg2(发电机2功率动作)功率变化100*xAg3(发电机3功率动作)功率变化200*xAv(太阳能电压动作):太阳能电压变为(1+x)倍。'
'三、功率计算公式'
'1.Pc(电池充放电功率)= max(0.2,min(0.8,Sc+0.18*Ac))2.Pg1(发电机1功率)= max(0,min(150,100*Ag1+Sg1))'
'3.Pg2(发电机2功率)= max(0,min(375,100*Ag2+Sg2))4.Pg3(发电机3功率)= max(0,min(500,200*Ag3+Sg3))'
'5.Pso(太阳能功率)= max(0,(0.45*Si+0.015*St60.3)*(1+Av))'
'6.Pw(风能功率) if 3<=Sw<8, Pw=Sw^3*172.2625/1000 elif 8<=Sw<12, Pw=64*172.2625/125 else Pw=0'
'四、成本和能源惩罚计算公式'
'P(总功率)=Pc+Pg1+Pg2+Pg3+Pso+Pw'
'if P>=Sl, Ee(过剩能源)= PSl, else Es(短缺能源)=Sl-P'
'Cb(电池成本)=0.01*Ac+0.1*Sc;Cg1(发电机1成本)=0.0034*Pg1^2+3*Pg1+30'
'Cg2(发电机2成本)=0.001*Pg2^2+10*Pg2+40;Cg3(发电机3成本)=0.001*Pg3^2+15*Pg3+70'
'Cs(太阳能成本)=0.01*Pso;Cw(风能成本)=0.01*Pw'
'Rs(销售收入)=0.5*Sp*Ee;Cp(购买成本)=Sp*Es'
'Pe(过剩惩罚)=(Ee-100)*50;Ps(短缺惩罚)=(Es-100)*50'
'Rs、Pe和Cp、Ps对应能源过剩和能源短缺情况不能同时存在。'
'五、目标函数 Goal=- (Cb+Cg1+Cg2+Cg3+Cs+Cw+Pe+Ps-Rs+Cp)/1000'
'六、控制方案和目标:'
'1.优先使用太阳能和风能发电,当无法满足负荷时,使用发电机和电池供电;'
'2.供电时综合考虑各发电模块的成本,通常前一台发电机的最大功率不足以满足负荷时,才开启后续发电机;'
'3.减少超过公网交易限额的交易,避免公共电网波动影响住宅用电稳定性;'
'4.满足下一时刻住宅用电负载,最大化奖励函数且使供需差异最小化(Ee、Es尽可能小),提升训练效果和策略质量。'
'请你读取每小时的Sp、Sl、St、Si和Sw数据并计算每个时刻的状态然后使用最合适的优化方法给出每小时的动作决策'
'其长度为5分别表示Ac、Ag1、Ag2、Ag3和Av每小时计算一次。'
'结果应为Json格式[{{x1, x2, x3, x4, x5}}, {{x1, x2, x3, x4, x5}}, ...],如果计算量太大,可以分批回答。'}],
temperature=0.8,
top_p=0.8
)
print(completion.model_dump_json())
if __name__ == '__main__':
get_response()

190
generate_optimization.py Normal file
View File

@ -0,0 +1,190 @@
import numpy as np
import pandas as pd
import json
from scipy.optimize import minimize
from pyswarm import pso
# 读取数据
file_path = './data.csv'
data = pd.read_csv(file_path)
# 初始化状态
initial_soc = 0.4 # 电池初始荷电状态
initial_sg1 = 0
initial_sg2 = 0
initial_sg3 = 0
battery_capacity = 500 * initial_soc
# 参数设置
battery_efficiency = 0.9
battery_degradation = 0.01
battery_holding = 0.1
solar_cofficient = 0.01
volatage_change_percent = 0.1
wind_cofficient = 0.01
# 发电机参数
DG_parameters = {
'gen_1': {'a': 0.0034, 'b': 3, 'c': 30, 'power_output_max': 150, 'power_output_min': 10, 'ramping_up': 100,
'ramping_down': 100},
'gen_2': {'a': 0.001, 'b': 10, 'c': 40, 'power_output_max': 375, 'power_output_min': 50, 'ramping_up': 100,
'ramping_down': 100},
'gen_3': {'a': 0.001, 'b': 15, 'c': 70, 'power_output_max': 500, 'power_output_min': 100, 'ramping_up': 200,
'ramping_down': 200}
}
NUM_GEN = len(DG_parameters.keys())
period = 100 # 处理前 100 个数据
def get_dg_info(parameters):
p_max = []
p_min = []
ramping_up = []
ramping_down = []
a_para = []
b_para = []
c_para = []
for name, gen_info in parameters.items():
p_max.append(gen_info['power_output_max'])
p_min.append(gen_info['power_output_min'])
ramping_up.append(gen_info['ramping_up'])
ramping_down.append(gen_info['ramping_down'])
a_para.append(gen_info['a'])
b_para.append(gen_info['b'])
c_para.append(gen_info['c'])
return p_max, p_min, ramping_up, ramping_down, a_para, b_para, c_para
p_max, p_min, ramping_up, ramping_down, a_para, b_para, c_para = get_dg_info(parameters=DG_parameters)
# 目标函数
def objective(x, period_batch, start, initial_states):
total_cost = 0
soc, sg1, sg2, sg3 = initial_states
for t in range(period_batch):
Ac, Ag1, Ag2, Ag3, Av = x[t * 5:(t + 1) * 5]
Pc = max(0.2, min(0.8, soc + 0.2 * Ac * 0.9))
Pg1 = max(0, min(150, 100 * Ag1 + sg1))
Pg2 = max(0, min(375, 100 * Ag2 + sg2))
Pg3 = max(0, min(500, 200 * Ag3 + sg3))
Pso = max(0, (0.2 * data['irradiance'][start + t] + 0.05 * data['temperature'][start + t] - 9.25) * (1 + Av))
if 3 <= data['wind_speed'][start + t] < 8:
Pw = data['wind_speed'][start + t] ** 3 * 172.2625 / 1000
elif 8 <= data['wind_speed'][start + t] < 12:
Pw = 64 * 172.2625 / 125
else:
Pw = 0
P = Pc + Pg1 + Pg2 + Pg3 + Pso + Pw
if P >= data['load'][start + t]:
Ee = P - data['load'][start + t]
Es = 0
else:
Es = data['load'][start + t] - P
Ee = 0
Cb = 0.01 * battery_capacity + 0.1 * Pc
Cg1 = 0.0034 * Pg1 ** 2 + 3 * Pg1 + 30
Cg2 = 0.001 * Pg2 ** 2 + 10 * Pg2 + 40
Cg3 = 0.001 * Pg3 ** 2 + 15 * Pg3 + 70
Cs = 0.01 * Pso
Cw = 0.01 * Pw
Rs = 0.5 * data['price'][start + t] * Ee
Cp = data['price'][start + t] * Es
Pe = (Ee - 100) * 50 if Ee > 100 else 0
Ps = (Es - 100) * 50 if Es > 100 else 0
total_cost += (Cb + Cg1 + Cg2 + Cg3 + Cs + Cw + Pe + Ps - Rs + Cp) / 1000
soc = Pc
sg1 = Pg1
sg2 = Pg2
sg3 = Pg3
return total_cost
# 约束条件
def constraint(x, period_batch, start, initial_states):
constraints = []
soc, sg1, sg2, sg3 = initial_states
for t in range(period_batch):
Ac, Ag1, Ag2, Ag3, Av = x[t * 5:(t + 1) * 5]
Pc = max(0.2, min(0.8, soc + 0.2 * Ac * 0.9))
Pg1 = max(10, min(150, 100 * Ag1 + sg1))
Pg2 = max(50, min(375, 100 * Ag2 + sg2))
Pg3 = max(100, min(500, 200 * Ag3 + sg3))
Pso = max(0, (0.2 * data['irradiance'][start + t] + 0.05 * data['temperature'][start + t] - 9.25) * (1 + Av))
if 3 <= data['wind_speed'][start + t] < 8:
Pw = data['wind_speed'][start + t] ** 3 * 172.2625 / 1000
elif 8 <= data['wind_speed'][start + t] < 12:
Pw = 64 * 172.2625 / 125
else:
Pw = 0
P = Pc + Pg1 + Pg2 + Pg3 + Pso + Pw
# 保证变化后的值不小于最小值
constraints.append(Pc - 0.2)
constraints.append(Pg1 - 10)
constraints.append(Pg2 - 50)
constraints.append(Pg3 - 100)
constraints.append(P - data['load'][start + t])
soc = Pc
sg1 = Pg1
sg2 = Pg2
sg3 = Pg3
return constraints
actions = []
batch_size = 10 # 每次处理 10 个数据
total_batches = 10 # 总共处理 100 个数据
start = 0
for batch_num in range(total_batches):
end = min(start + batch_size, period)
period_batch = end - start
# 初始化初始状态
initial_states = (initial_soc, initial_sg1, initial_sg2, initial_sg3)
# # 使用随机初始化更接近实际情况
# x0 = np.random.uniform(0, 1, period_batch * 5)
# bounds = [(-1, 1)] * (period_batch * 5)
# cons = {'type': 'ineq', 'fun': lambda x: constraint(x, period_batch, start, initial_states)}
# result = minimize(lambda x: objective(x, period_batch, start, initial_states), x0, method='SLSQP', bounds=bounds,
# constraints=cons)
#
# actions.extend(result.x.reshape((period_batch, 5)).tolist())
# 使用粒子群优化算法
lb = [-1] * (period_batch * 5)
ub = [1] * (period_batch * 5)
xopt, fopt = pso(lambda x: objective(x, period_batch, start, initial_states), lb, ub,
ieqcons=[lambda x: constraint(x, period_batch, start, initial_states)], maxiter=100, swarmsize=50)
actions.extend(xopt.reshape((period_batch, 5)).tolist())
# 更新初始状态为当前批次最后一个状态
last_action = actions[-1]
initial_soc = max(0.2, min(0.8, initial_soc + 0.2 * last_action[0] * 0.9))
initial_sg1 = max(10, min(150, 100 * last_action[1] + initial_sg1))
initial_sg2 = max(50, min(375, 100 * last_action[2] + initial_sg2))
initial_sg3 = max(100, min(500, 200 * last_action[3] + initial_sg3))
start += batch_size
actions_json_path = './actions_batch_p.json'
try:
with open(actions_json_path, 'w') as f:
json.dump(actions, f)
print(f"文件保存成功:{actions_json_path}")
except Exception as e:
print(f"文件保存失败:{e}")

183
genetic_gpu_all.py Normal file
View File

@ -0,0 +1,183 @@
import json
import os
import time
import numpy as np
import pandas as pd
import torch
from concurrent.futures import ThreadPoolExecutor
from deap import base, creator, tools, algorithms
def fitness_torch(individuals, price, load, temperature, irradiance, wind_speed, device):
individuals_torch = torch.tensor(individuals, device=device)
num_individuals = individuals_torch.shape[0]
Ac, Ag1, Ag2, Ag3, Av = [individuals_torch[:, i * period:(i + 1) * period] for i in range(5)]
soc = torch.zeros((num_individuals, period), device=device)
Pg1, Pg2, Pg3 = [torch.zeros((num_individuals, period), device=device) for _ in range(3)]
Cb, Cg1, Cg2, Cg3, Cs, Cw, Rs, Cp, Pe, Ps, Ee, Es = [torch.zeros((num_individuals, period), device=device) for _ in
range(12)]
price_torch = torch.tensor(price, device=device)
load_torch = torch.tensor(load, device=device)
temperature_torch = torch.tensor(temperature, device=device)
irradiance_torch = torch.tensor(irradiance, device=device)
wind_speed_torch = torch.tensor(wind_speed, device=device)
for t in range(period):
if t > 0:
soc[:, t] = torch.clamp(soc[:, t - 1] + 0.2 * Ac[:, t] * 0.9, 0.2, 0.8)
Pg1[:, t] = torch.clamp(Pg1[:, t - 1] + 100 * Ag1[:, t], min=0, max=150)
Pg2[:, t] = torch.clamp(Pg2[:, t - 1] + 100 * Ag2[:, t], min=0, max=375)
Pg3[:, t] = torch.clamp(Pg3[:, t - 1] + 200 * Ag3[:, t], min=0, max=500)
else:
soc[:, t] = 0.4
Pg1[:, t] = torch.clamp(100 * Ag1[:, t], min=0, max=150)
Pg2[:, t] = torch.clamp(100 * Ag2[:, t], min=0, max=375)
Pg3[:, t] = torch.clamp(200 * Ag3[:, t], min=0, max=500)
Pso = torch.clamp((0.2 * irradiance_torch[t] + 0.05 * temperature_torch[t] - 9.25) * (1 + Av[:, t]), min=0)
Pw = torch.where(
(wind_speed_torch[t] >= 3) & (wind_speed_torch[t] < 8),
wind_speed_torch[t] ** 3 * 172.2625 / 1000,
torch.where(
(wind_speed_torch[t] >= 8) & (wind_speed_torch[t] < 12),
64 * 172.2625 / 125,
torch.zeros_like(wind_speed_torch[t])
)
)
P = Ac[:, t] + Pg1[:, t] + Pg2[:, t] + Pg3[:, t] + Pso + Pw
Ee[:, t] = torch.where(P >= load_torch[t], P - load_torch[t], torch.zeros_like(P))
Es[:, t] = torch.where(P < load_torch[t], load_torch[t] - P, torch.zeros_like(P))
Cb[:, t] = 0.01 * Ac[:, t] + 0.1 * soc[:, t]
Cg1[:, t] = 0.0034 * Pg1[:, t] ** 2 + 3 * Pg1[:, t] + 30
Cg2[:, t] = 0.001 * Pg2[:, t] ** 2 + 10 * Pg2[:, t] + 40
Cg3[:, t] = 0.001 * Pg3[:, t] ** 2 + 15 * Pg3[:, t] + 70
Cs[:, t] = 0.01 * Pso
Cw[:, t] = 0.01 * Pw
Rs[:, t] = 0.5 * price_torch[t] * Ee[:, t]
Cp[:, t] = price_torch[t] * Es[:, t]
Pe[:, t] = torch.where(Ee[:, t] > 100, (Ee[:, t] - 100) * 50, torch.zeros_like(Ee[:, t]))
Ps[:, t] = torch.where(Es[:, t] > 100, (Es[:, t] - 100) * 50, torch.zeros_like(Es[:, t]))
total_cost = torch.sum(Cb + Cg1 + Cg2 + Cg3 + Cs + Cw + Pe + Ps - Rs + Cp, dim=1)
reward = -total_cost / 1000
return reward.cpu().numpy()
def check_bounds(func):
def wrapper(*args, **kwargs):
offspring = func(*args, **kwargs)
for individual in offspring:
for i in range(len(individual)):
individual[i] = np.clip(individual[i], -1, 1)
return offspring
return wrapper
def save_progress(population, gen, filename='ga_progress.json'):
data = {
'population': [[list(ind), ind.fitness.values[0]] for ind in population],
'generation': gen
}
with open(filename, 'w') as f:
json.dump(data, f)
print(f"进度已保存到第 {gen} 代到 '{filename}'")
def load_progress(filename='ga_progress.json'):
if os.path.exists(filename):
with open(filename, 'r') as f:
data = json.load(f)
population = []
for ind_data in data['population']:
ind = creator.Individual(ind_data[0])
ind.fitness.values = (ind_data[1],)
population.append(ind)
gen = data['generation']
print(f"从第 {gen} 代加载进度。")
return population, gen
else:
return None, 0
def save_decision_values(best_ind, period, file_suffix):
decision_values = [
{"x1": best_ind[i], "x2": best_ind[i + period], "x3": best_ind[i + 2 * period], "x4": best_ind[i + 3 * period],
"x5": best_ind[i + 4 * period]} for i in range(period)]
filename = f'./decision_values_{file_suffix}.json'
with open(filename, 'w') as f:
json.dump(decision_values, f)
print(f"周期 {file_suffix}:决策值已保存到 '{filename}'")
def main():
data = pd.read_csv('./data.csv')
# period = len(data)
period = 2
price, load, temperature, irradiance, wind_speed = [data[col].values for col in
['price', 'load', 'temperature', 'irradiance', 'wind_speed']]
creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("Individual", list, fitness=creator.FitnessMax)
toolbox = base.Toolbox()
toolbox.register("attr_float", np.random.uniform, -1, 1)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, n=5 * period)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("mate", check_bounds(tools.cxBlend), alpha=0.5)
toolbox.register("mutate", check_bounds(tools.mutGaussian), mu=0, sigma=0.1, indpb=0.2)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("evaluate", fitness_torch)
population, start_gen = load_progress()
if population is None:
population = toolbox.population(n=500)
start_gen = 0
NGEN, CXPB, MUTPB = 500, 0.7, 0.2
batch_size = 500
for gen in range(start_gen, NGEN):
start_time = time.time()
offspring = algorithms.varAnd(population, toolbox, cxpb=CXPB, mutpb=MUTPB)
num_individuals = len(offspring)
with ThreadPoolExecutor() as executor:
futures = []
for i in range(0, num_individuals, batch_size):
batch = offspring[i:i + batch_size]
individuals = [ind[:] for ind in batch]
futures.append(
executor.submit(toolbox.evaluate, individuals, price, load, temperature, irradiance, wind_speed, 0))
for future in futures:
fitnesses = future.result()
for ind, fitness in zip(offspring, fitnesses):
ind.fitness.values = (fitness,)
population = toolbox.select(offspring, k=len(population))
end_time = time.time()
elapsed_time = end_time - start_time
print(f"{gen + 1} 代完成于 {elapsed_time:.2f}")
if (gen + 1) % 100 == 0:
best_ind = tools.selBest(population, 1)[0]
save_decision_values(best_ind, period, gen + 1)
save_progress(population, gen + 1)
best_ind = tools.selBest(population, 1)[0]
print('最佳个体:', best_ind)
print('适应度:', best_ind.fitness.values)
save_decision_values(best_ind, period, NGEN)
save_progress(population, NGEN)
if __name__ == "__main__":
main()

183
genrtic_gpu_7.22.py Normal file
View File

@ -0,0 +1,183 @@
import json
import os
import time
import numpy as np
import pandas as pd
import torch
from concurrent.futures import ThreadPoolExecutor
from deap import base, creator, tools, algorithms
def fitness_torch(individuals, price, load, temperature, irradiance, wind_speed, prev_soc, prev_Pg1, prev_Pg2, prev_Pg3,
device):
individuals_torch = torch.tensor(individuals, device=device)
num = individuals_torch.shape[0]
Ac, Ag1, Ag2, Ag3, Av = [individuals_torch[:, i * period:(i + 1) * period] for i in range(5)]
# soc = torch.zeros((num, period), device=device)
# Pg1, Pg2, Pg3 = [torch.zeros((num, period), device=device) for _ in range(3)]
# Cb, Cg1, Cg2, Cg3, Cs, Cw, Rs, Cp, Pe, Ps, Ee, Es = [torch.zeros((num, period), device=device) for _ in range(12)]
soc = torch.clamp(prev_soc + 0.2 * Ac * 0.9, 0.2, 0.8)
Pg1 = torch.clamp(prev_Pg1 + 100 * Ag1, min=0, max=150)
Pg2 = torch.clamp(prev_Pg2 + 100 * Ag2, min=0, max=375)
Pg3 = torch.clamp(prev_Pg3 + 200 * Ag3, min=0, max=500)
Pso = torch.clamp((0.2 * irradiance + 0.05 * temperature - 9.25) * (1 + Av), min=0)
Pw = torch.where(
(wind_speed >= 3) & (wind_speed < 8),
wind_speed ** 3 * 172.2625 / 1000,
torch.where(
(wind_speed >= 8) & (wind_speed < 12),
64 * 172.2625 / 125,
torch.zeros_like(wind_speed)
)
)
P = Ac + Pg1 + Pg2 + Pg3 + Pso + Pw
Ee = torch.where(P >= load, P - load, torch.zeros_like(P))
Es = torch.where(P < load, load - P, torch.zeros_like(P))
Cb = 0.01 * Ac + 0.1 * soc
Cg1 = 0.0034 * Pg1 ** 2 + 3 * Pg1 + 30
Cg2 = 0.001 * Pg2 ** 2 + 10 * Pg2 + 40
Cg3 = 0.001 * Pg3 ** 2 + 15 * Pg3 + 70
Cs = 0.01 * Pso
Cw = 0.01 * Pw
Rs = 0.5 * price * Ee
Cp = price * Es
Pe = torch.where(Ee > 100, (Ee - 100) * 50, torch.zeros_like(Ee))
Ps = torch.where(Es > 100, (Es - 100) * 50, torch.zeros_like(Es))
total_cost = torch.sum(Cb + Cg1 + Cg2 + Cg3 + Cs + Cw + Pe + Ps - Rs + Cp, dim=1)
reward = -total_cost / 1000
return reward.cpu().numpy(), soc.cpu().numpy(), Pg1.cpu().numpy(), Pg2.cpu().numpy(), Pg3.cpu().numpy()
def save_decision_values(best_ind, period, index):
decisions = {
'Ac': best_ind[0],
'Ag1': best_ind[1],
'Ag2': best_ind[2],
'Ag3': best_ind[3],
'Av': best_ind[4]
}
with open(f'decision_values_{period}_index_{index}.json', 'w') as f:
json.dump(decisions, f)
def save_progress(population, period):
population_data = [ind.tolist() for ind in population]
with open(f'population_{period}.json', 'w') as f:
json.dump({'period': period, 'population': population_data}, f)
def load_progress():
if os.path.exists('population_gen_499.json'):
with open('population_gen_499.json', 'r') as f:
data = json.load(f)
return data['population'], data['period']
return None, 0
def check_bounds(func):
def wrapper(*args, **kwargs):
offspring = func(*args, **kwargs)
if offspring[0] is None or offspring[1] is None:
print("Error: One of the offspring is None", offspring)
raise ValueError("Offspring cannot be None.")
for child in offspring:
for i in range(len(child)):
if child[i] < -1:
child[i] = -1
elif child[i] > 1:
child[i] = 1
return offspring
return wrapper
def main():
period = 8760
NGEN = 500
CXPB, MUTPB = 0.7, 0.2
batch_size = 500
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
data = pd.read_csv('./data.csv')
price, load, temperature, irradiance, wind_speed = [data[col].values for col in
['price', 'load', 'temperature', 'irradiance', 'wind_speed']]
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", list, fitness=creator.FitnessMin)
toolbox = base.Toolbox()
toolbox.register("attr_float", np.random.uniform, -1, 1)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, 5)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("mate", check_bounds(tools.cxBlend), alpha=0.5)
toolbox.register("mutate", check_bounds(tools.mutGaussian), mu=0, sigma=0.1, indpb=0.2)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("evaluate", fitness_torch)
population = load_progress()
if population is None:
population = toolbox.population(n=NGEN)
print("Initial population:", len(population), "and first individual:", population[0])
prev_soc = torch.tensor([0.4] * NGEN, device=device)
prev_Pg1 = torch.zeros(NGEN, device=device)
prev_Pg2 = torch.zeros(NGEN, device=device)
prev_Pg3 = torch.zeros(NGEN, device=device)
for index in range(period):
for gen in range(NGEN):
start_time = time.time()
offspring = algorithms.varAnd(population, toolbox, cxpb=CXPB, mutpb=MUTPB)
num_individuals = len(offspring)
with ThreadPoolExecutor() as executor:
futures = []
for i in range(0, num_individuals, batch_size):
batch = offspring[i:i + batch_size]
individuals = [ind[:] for ind in batch]
futures.append(
executor.submit(toolbox.evaluate, individuals, price[index], load[index], temperature[index],
irradiance[index], wind_speed[index], prev_soc, prev_Pg1, prev_Pg2, prev_Pg3,
device)
)
for future in futures:
fitnesses, socs, Pg1s, Pg2s, Pg3s = zip(*future.result())
for ind, fitness in zip(offspring, fitnesses):
ind.fitness.values = (fitness,)
prev_soc[:len(socs)] = torch.tensor(socs, device=device)
prev_Pg1[:len(Pg1s)] = torch.tensor(Pg1s, device=device)
prev_Pg2[:len(Pg2s)] = torch.tensor(Pg2s, device=device)
prev_Pg3[:len(Pg3s)] = torch.tensor(Pg3s, device=device)
population = toolbox.select(offspring, k=len(population))
print("Population after selection:", population[:5])
end_time = time.time()
print(f"{index + 1}小时完成花费 {end_time - start_time:.2f}")
best_ind = tools.selBest(population, 1)[0]
print('最佳个体:', best_ind)
print('适应度:', best_ind.fitness.values)
save_decision_values(best_ind, period, index)
save_progress(population, period)
prev_soc.fill_(0.4)
prev_Pg1.fill_(0)
prev_Pg2.fill_(0)
prev_Pg3.fill_(0)
prev_soc[:len(best_ind)] = torch.tensor(best_ind[:period], device=device)
prev_Pg1[:len(best_ind)] = torch.tensor(best_ind[period:2 * period], device=device)
prev_Pg2[:len(best_ind)] = torch.tensor(best_ind[2 * period:3 * period], device=device)
prev_Pg3[:len(best_ind)] = torch.tensor(best_ind[3 * period:4 * period], device=device)
if __name__ == "__main__":
main()

50
llm.txt Normal file
View File

@ -0,0 +1,50 @@
目前的能源系统结构主要由太阳能发电站、风能发电站、发电机和电池模块通过电力网与住宅相连接。以下是中文和符号的对照:
As电池的充电状态动作Ag1发电机1的功率动作Ag2发电机2的功率动作Ag3发电机3的功率动作Av太阳能的电压动作
Sp电价状态Sc电池荷电状态Sl住宅负载状态Sg1发电机1功率状态Sg2发电机2功率状态Sg3发电机3功率状态St温度状态Si太阳辐射度状态Sw风速状态
Pc电池充放电功率Pg1发电机1功率Pg2发电机2功率Pg3发电机3功率Pso太阳能功率Pw风能功率
Cb电池成本Cg1发电机1成本Cg2发电机2成本Cg3发电机3成本Cs太阳能成本Cw风能成本Rs销售收入Cp购买成本Ee过剩能源Es短缺能源Pe过剩惩罚Ps短缺惩罚。
一、控制方案和目标:
1、优先使用太阳能和风能发电当无法满足满负荷时请使用发电机和电池供电。
2、在供电时考虑多种供电来源太阳能、风能、发电机和电池综合每个发电模块的成本需满足下一时刻的住宅用电负载且使供需差异最小化。
3、减少超过公网交易限额的交易以避免公共电网波动影响住宅用电的稳定性。
4、最大化奖励函数并最小化能源不平衡。
二、连续状态空间是一个长度为9的向量包含以下部分
Sp由数据提供Sc初始为0.4始终介于0.2和0.8Sl由数据提供
Sg1初始功率为0Sg2初始功率为0Sg3初始功率为0
St由数据提供Si由数据提供Sw由数据提供。
其中ScSg1、Sg2和Sg3始终不小于0。
三、连续动作空间一个长度为5的向量x表示当前时间的取值, 其范围从-1到1每个时刻的值都是在上一时刻基础上进行变化的Ag1、Ag2和Ag3分别表示相对于各自最大瞬时功率100、100和200的变化值。
1、Ac电池总容量为500初始电池容量为500*0.4=200电池容量变化100*x
2、Ag1发电量变化100*x
3、Ag2发电量变化100*x
4、Ag3发电量变化200*x
5、Av太阳能电压变为(1+x)倍。
可根据需求通过改变功率控制发电机的开关当关闭时保持Sg=0当打开时需满足Pg约束。根据成本通常前一台发电机的最大功率不足以满足剩余负荷时才开启后续发电机。
四、发电量计算公式
1、Pc = max(0.2, min(0.8, Sc + 0.2 * action[0] * 0.9))其中0.8和0.2为电池荷电状态上下限;
2、Pg1 = max(0, min(150, 100 * action[1] + Sg1))其中150和0为发电机1开启时发电量上下限
3、Pg2 = max(0, min(375, 100 * action[2] + Sg2))其中375和0为发电机2开启时发电量上下限
4、Pg3 = max(0, min(500, 200 * action[3] + Sg3))其中500和0为发电机2开启时发电量上下限
5、Pso = max(0, (0.2 * Si + 0.05 * St - 9.25) * (1 + action[4]))其中太阳能发电功率不小于0。
6、Pw与Ws有关
if 3 <= Ws < 8, Pw = Ws^3 * 172.2625 / 1000
elif 8 <= Ws < 12, Pw = 64 * 172.2625 / 125
else Pw = 0
五、成本和能源平衡计算
P = Pc + Pg1 + Pg2 + Pg3 + Pso + Pw
if P >= Sl , Ee = P Sl, else Es = Sl - P
1、Cb = 0.01 * Asoc + 0.1 * Sc
2、Cg1 = 0.0034 * Pg1^2 + 3 * Pg1 + 30
3、Cg2 = 0.001 * Pg2^2 + 10 * Pg2 + 40
4、Cg3 = 0.001 * Pg3^2 + 15 * Pg3 + 70
5、Cs = 0.01 * Pso
6、Cw = 0.01 * Pw
7、Rs = 0.5 * Sp * Ee
8、Cp = Sp * Es
9、Pe = (Ee - 100) * 50
10、Ps = (Es - 100) * 50
因为Rs、Pe和Cp、Ps对应能源过剩和能源短缺情况只能同时存在一组。
六、奖励函数
Reward = - (Cb + Cg1 + Cg2 + Cg3 + Cs + Cw + Pe + Ps - Rs + Cp) / 1000
七、请求动作数据
读取数据并使用上述公式计算每个时刻的状态请你使用最合适的优化方法给出每个时刻的动作控制方案每小时计算一次并生成一个长度为5的数组表示Ac、Ag1、Ag2、Ag3和Av的动作决策。返回的数据应为以下Json格式[{x1, x2, x3, x4, x5}, {x1, x2, x3, x4, x5}, ..., {x1, x2, x3, x4, x5}]

View File

@ -39,7 +39,7 @@ llm = ChatOpenAI(
verbose=True,
openai_api_key="none",
# openai_api_base="http://0.0.0.0:5049/v1/models",
openai_api_base="http://localhost:8000/v1",
openai_api_base="http://0.0.0.0:8501",
model_name="Qwen1.5-32b-int4"
)
prompt = ChatPromptTemplate.from_messages(
@ -65,7 +65,7 @@ for i in range(num_hours):
csv_data_chunk = format_csv_data(df, start, end)
result = with_message_history.invoke(
{"input": f"数据如下:\n{csv_data_chunk}\n只返回json格式的五个决策数据{{[x1 x2 x3 x4 x5]}}"},
{"input": f"数据如下:\n{csv_data_chunk}"},
config={
"configurable": {"session_id": "cxd"}
},

View File

@ -58,7 +58,7 @@ class Battery:
self.current_capacity = updated_capacity # update capacity to current state
def get_cost(self, energy_change, energy_hold): # cost depends on the energy change
cost = energy_change * self.degradation + energy_hold * self.holding
cost = abs(energy_change) * self.degradation + energy_hold * self.holding
return cost
def SOC(self):

View File

@ -1,5 +1,5 @@
solar_parameters = {
'I_sc0': 8.0, # 参考条件下的短路电流 (A)
'I_sc0': 10.0, # 参考条件下的短路电流 (A)
'V_b': 25, # 基准电压
'V_oc0': 36.0, # 参考条件下的开路电压 (V)
'R_s': 0.1, # 串联电阻 (Ω)

88
qwen_connect.py Normal file
View File

@ -0,0 +1,88 @@
from pathlib import Path
from openai import OpenAI
import os
client = OpenAI(
api_key='sk-f01744b2801344b1a72f89ec7e290cad',
base_url="https://dashscope.aliyuncs.com/compatible-mode/v1", # 填写DashScope服务base_url
)
file_csv = client.files.create(file=Path("./data.csv"), purpose="file-extract")
# 新文件上传后需要等待模型解析,首轮响应时间可能较长
completion = client.chat.completions.create(
model="qwen-long",
messages=[
{
'role': 'system',
'content':
'你是一名擅长调试和控制智慧住宅能源系统的大模型和强化学习专家,可根据当前的环境状态给出最优的动作决策。'
'目前的能源系统结构主要由太阳能发电站、风能发电站、发电机和电池模块通过电力网与住宅相连接。以下是中文和符号的对照:'
'As电池的充电状态动作Ag1发电机1的功率动作Ag2发电机2的功率动作Ag3发电机3的功率动作Av太阳能的电压动作'
'Sp电价状态Sc电池荷电状态Sl住宅负载状态Sg1发电机1功率状态Sg2发电机2功率状态Sg3发电机3功率状态St温度状态Si太阳辐射度状态Sw风速状态'
'Pc电池充放电量Pg1发电机1发电量Pg2发电机2发电量Pg3发电机3发电量Pso太阳能发电量Pw风能发电量'
'Cb电池成本Cg1发电机1成本Cg2发电机2成本Cg3发电机3成本Cs太阳能成本Cw风能成本'
'Rs销售收入Cp购买成本Ee过剩能源Es短缺能源Pe过剩惩罚Ps短缺惩罚。'
'一、控制方案和目标:'
'1、优先使用太阳能和风能发电当无法满足满负荷时请使用发电机和电池供电。'
'2、在供电时考虑多种供电来源太阳能、风能、发电机和电池综合每个发电模块的成本需满足下一时刻的住宅用电负载且使供需差异最小化。'
'3、减少超过公网交易限额的交易以避免公共电网波动影响住宅用电的稳定性。'
'4、最大化奖励函数并最小化能源不平衡。'
'二、连续状态空间是一个长度为9的向量包含以下部分'
'Sp由数据提供Sc初始为0.4始终介于0.2和0.8Sl由数据提供'
'Sg1初始功率为0Sg2初始功率为0Sg3初始功率为0'
'St由数据提供Si由数据提供Sw由数据提供。'
'其中ScSg1、Sg2和Sg3始终不小于0。'
'三、连续动作空间一个长度为5的向量x表示当前时间的取值, 其范围从-1到1'
'每个时刻的值都是在上一时刻基础上进行变化的Ag1、Ag2和Ag3分别表示相对于各自最大瞬时功率100、100和200的变化值。'
'1、Ac电池总容量为500初始电池容量为500*0.4=200电池容量变化100*x'
'2、Ag1发电量变化100*x'
'3、Ag2发电量变化100*x'
'4、Ag3发电量变化200*x'
'5、Av太阳能电压变为(1+x)倍。'
'可根据需求通过改变功率控制发电机的开关当关闭时保持Sg=0当打开时需满足Pg约束。'
'根据成本,通常前一台发电机的最大功率不足以满足剩余负荷时,才开启后续发电机。'
'四、发电量计算公式'
'1、Pc = max(0.2, min(0.8, Sc + 0.2 * action[0] * 0.9))其中0.8和0.2为电池荷电状态上下限;'
'2、Pg1 = max(0, min(150, 100 * action[1] + Sg1))其中150和0为发电机1开启时发电量上下限'
'3、Pg2 = max(0, min(375, 100 * action[2] + Sg2))其中375和0为发电机2开启时发电量上下限'
'4、Pg3 = max(0, min(500, 200 * action[3] + Sg3))其中500和0为发电机2开启时发电量上下限'
'5、Pso = max(0, (0.2 * Si + 0.05 * St - 9.25) * (1 + action[4]))其中太阳能发电功率不小于0。'
'6、Pw与Ws有关'
'if 3 <= Ws < 8, Pw = Ws^3 * 172.2625 / 1000'
'elif 8 <= Ws < 12, Pw = 64 * 172.2625 / 125'
'else Pw = 0'
'五、成本和能源平衡计算'
'P = Pc + Pg1 + Pg2 + Pg3 + Pso + Pw'
'if P >= Sl , Ee = P Sl, else Es = Sl - P'
'1、Cb = 0.01 * Asoc + 0.1 * Sc'
'2、Cg1 = 0.0034 * Pg1^2 + 3 * Pg1 + 30'
'3、Cg2 = 0.001 * Pg2^2 + 10 * Pg2 + 40'
'4、Cg3 = 0.001 * Pg3^2 + 15 * Pg3 + 70'
'5、Cs = 0.01 * Pso'
'6、Cw = 0.01 * Pw'
'7、Rs = 0.5 * Sp * Ee'
'8、Cp = Sp * Es'
'9、Pe = (Ee - 100) * 50'
'10、Ps = (Es - 100) * 50'
'因为Rs、Pe和Cp、Ps对应能源过剩和能源短缺情况只能同时存在一组。'
'六、奖励函数'
'Reward = - (Cb + Cg1 + Cg2 + Cg3 + Cs + Cw + Pe + Ps - Rs + Cp) / 1000'
},
{
'role': 'system',
'content': f'fileid://{file_csv.id}'
},
{
'role': 'user',
'content':
'读取数据并使用上述公式计算每个时刻的状态,请你使用最合适的优化方法,给出每个时刻的动作控制方案,'
'每小时计算一次并生成一个长度为5的数组表示Ac、Ag1、Ag2、Ag3和Av的动作决策。'
'返回的数据应为以下Json格式[{x1, x2, x3, x4, x5}, {x1, x2, x3, x4, x5}, ..., {x1, x2, x3, x4, x5}]'
'只需要json数据不需要其他任何分析。'
}
],
stream=False
)
print(completion.choices[0].message.model_dump()) # 非流式输出

18
test.py
View File

@ -52,9 +52,17 @@
#
# data = get_llm_action(2)
# print(data)
from environment import ESSEnv
import torch
env = ESSEnv()
wind_speed = env.data_manager.get_series_wind_data(1, 1)
wind = [env.wind.step(ws) for ws in wind_speed]
print(wind)
def get_available_gpus():
if torch.cuda.is_available():
num_gpus = torch.cuda.device_count()
print(f"Number of available GPUs: {num_gpus}")
for i in range(num_gpus):
print(f"GPU {i}: {torch.cuda.get_device_name(i)}")
else:
print("No GPUs are available.")
get_available_gpus()

View File

@ -41,7 +41,6 @@ def optimization_base_result(env, month, day, initial_soc):
battery_capacity = env.battery.capacity
battery_efficiency = env.battery.efficiency
solar_cofficient = env.solar.opex_cofficient
volatage_change_percent = env.solar.change_percent
wind_cofficient = env.wind.opex_cofficient
battery_degradation = env.battery.degradation
battery_holding = env.battery.holding
@ -62,8 +61,7 @@ def optimization_base_result(env, month, day, initial_soc):
pv_voltage = m.addVars(period, vtype=GRB.CONTINUOUS, lb=-1, ub=1, name='pv_voltage')
# 计算光伏和风力发电量
pv = [(0.2 * irradiance[t] + 0.05 * temperature[t] - 9.25) *
(1 + volatage_change_percent * pv_voltage[t]) for t in range(period)]
pv = [(0.2 * irradiance[t] + 0.05 * temperature[t] - 9.25) * (1 + pv_voltage[t]) for t in range(period)]
wind = [172.265625 * wind_speed[t] ** 3 / 1e3 if 3 <= wind_speed[t] < 8
else (172.265625 * 8 ** 3 / 1e3 if 8 <= wind_speed[t] < 12 else 0) for t in range(period)]
@ -99,8 +97,7 @@ def optimization_base_result(env, month, day, initial_soc):
m.optimize()
# 记录数据便于绘图
pv = [(0.2 * irradiance[t] + 0.05 * temperature[t] - 9.25) *
(1 + volatage_change_percent * pv_voltage[t].x) for t in range(period)]
pv = [(0.2 * irradiance[t] + 0.05 * temperature[t] - 9.25) * (1 + pv_voltage[t].x) for t in range(period)]
output_record = {'pv': [], 'wind': [], 'price': [], 'load': [], 'netload': [],
'soc': [], 'battery_energy_change': [], 'grid_import': [], 'grid_export': [],
'gen1': [], 'gen2': [], 'gen3': [], 'step_cost': []}