YuXT/01_Swarm_intelligence_optimization_algorithm
0
1import gradio as gr2import numpy as np3import pandas as pd4import math5from matplotlib import pyplot as plt6from math import ceil7 8 9# 读取文件10s = 111T_max = 50012count=48013 14tzz=pd.read_excel('data.xlsx', header=None)15Load = tzz.iloc[:count, 0] # 给定负荷16bat_power = tzz.iloc[0:count, 1]17pv_power = [0] * count18C_buy = tzz.iloc[0:count, 2] # 贵分时电价19 20 21# 如果上传文件，则读取22def choose_file(default_file, uploaded_file):23 global Load, bat_power, pv_power, C_buy24 if default_file:25 result, pic_path = main_func(20)26 27 else:28 filename = uploaded_file.name29 tzz = pd.read_excel(filename, header=None)30 Load = tzz.iloc[:count, 0] # 给定负荷31 bat_power = tzz.iloc[0:count, 1]32 pv_power = [0] * count33 C_buy = tzz.iloc[0:count, 2] # 贵分时电价34 result, pic_path = main_func(20)35 return result, pic_path36 37 38 39# 计算适应度函数的函数40def fitness_caculate(pp):41 # 参数42 duration = 0.0543 net_power = np.zeros(count)44 bat_cha_power = np.zeros(count)45 bat_dis_power = np.zeros(count)46 dim = 480 # dim表示粒子维度47 # 参数结束48 t_duration = 0.0549 F1 = 0 #F是惩罚项50 F2 = 051 F3 = 052 C4 = 0 #购电成本53 bat = pp[dim] * 0.454 bat_0 = pp[dim]55 bat_power = pp[0:count]56 T_max = 50057 rate = 0.0358 year = 2059 n = 22.68 * T_max / 3060 r = (rate * math.pow(1 + rate , year)) / (math.pow(1 + rate , year) - 1)61 sff = rate / (math.pow(1 + rate , year) - 1)62 sff_bat1 = rate / (math.pow(1 + rate , 6) - 1)63 sff_bat2 = rate / (math.pow(1 + rate , 12) - 1)64 sff_bat3 = rate / (math.pow(1 + rate , 18) - 1)65 sff_dc1 = rate / (math.pow(1 + rate , 8) - 1)66 sff_dc2 = rate / (math.pow(1 + rate , 16) - 1)67 sff_ac = rate / (math.pow(1 + rate , 10) - 1)68 S_bat = 2 / 369 S_dc = 1 / 270 S_ac = 071 for j in range(480):72 net_power[j] = (Load[j] + bat_power[j])73 if bat_power[j] > 0:74 bat_cha_power[j] = bat_power[j]75 bat_dis_power[j] = 076 bat = (bat + bat_cha_power[j] * t_duration)77 elif bat_power[j] == 0:78 bat_cha_power[j] = 079 bat_dis_power[j] = 080 bat = (bat + bat_dis_power[j] * t_duration)81 else:82 bat_cha_power[j] = 083 bat_dis_power[j] = bat_power[j]84 bat = (bat + bat_dis_power[j] * t_duration)85 86 if (0.3 * bat_0 <= bat) and (bat <= 0.9 * bat_0):87 F1 = F1 + 088 else:89 F1 = F1 + 10000 # F1用来惩罚储能容量超出l.90 91 if (net_power[j] <= T_max):92 F3 = F3 + 093 else:94 F3 = F3 + 1000095 C4 += net_power[j]* C_buy[j] * t_duration96 97 if abs(sum((1 / 0.95) * bat_dis_power[1: j]*t_duration + 0.95 * bat_cha_power[1: j]*t_duration)) <= 10:98 F2 = F2+099 else:100 F2 = F2+10000101 102 dc = max(abs(bat_power))103 dc = ceil(dc)104 ac = max(abs(net_power))105 ac = ceil(ac)106 C1 = r * (4800 * T_max + 1414.2 * bat_0 + 696.47 * ac + 471.56 * dc) / 365107 C2 = 0.01 * C1108 C3 = (1414.2 * bat_0 * (sff_bat1 + sff_bat2 + sff_bat3 - S_bat * sff) + \109 471.56 * dc * (sff_dc1+sff_dc2-S_dc * sff) + 696.47 * ac * (sff_ac - S_ac * sff)) / 365110 C4111 C = C4+ C1 + C2 + C3+ F1 + F2 + F3 + n112 113 return C4+ C1 + C2 + C3+ F1 + F2 + F3 + n114 115 116 117 118def main_func(sizepop):119 # 参数120 w_max = 0.9 # w越大，全局搜索能力越强，w越小，局部搜索能力越佳121 w_min = 0.4 # 我们用的权重是变化的122 v1_max = 10 # 粒子最大速度，设置速度限制123 v1_min = -10 # %如果粒子飞行速度过快，很可能直接飞过最优解位置，124 v2_max = 50 # 粒子最大速度，设置速度限制125 v2_min = -50 # %如果粒子飞行速度过快，很可能直接飞过最优解位置，126 # 设置种群参数127 sizepop = 20 # 粒子个数,每个粒子都是一个解128 iteration = 100 # 迭代次数129 dim = 480 # dim表示粒子维度130 B_max = 1300131 B_min = 1000132 # 参数结束133 # print('sizepop', sizepop)134 sizepop = int(sizepop)135 # 所有粒子位置和速度初始化136 pop_v = np.zeros((sizepop, dim + 1))137 pop_x = np.zeros((sizepop, dim + 1))138 for i in range(sizepop):139 pop_v[i, count] = v1_min + (v1_max - v1_min) * np.random.rand()140 pop_x[i, count] = B_min + (B_max - B_min) * np.random.rand()141 pop_x[i, count] = np.round(pop_x[i, count])142 for j in range(count):143 pop_v[i, j] = v2_min + (v2_max - v2_min) * np.random.rand()144 pop_x[i, j] = bat_power[j]145 146 # 保存最佳粒子的位置147 pbest_all = np.zeros((sizepop, dim + 1))148 pbest = np.zeros((sizepop, 1))149 150 151 # 参数结尾152 for i in range(sizepop):153 pop_v[i, count] = v1_min + (v1_max - v1_min) * np.random.rand()154 pop_x[i, count] = B_min + (B_max - B_min) * np.random.rand()155 pop_x[i, count] = np.round(pop_x[i, count])156 for j in range(count):157 pop_v[i, j] = v2_min + (v2_max - v2_min) * np.random.rand()158 pop_x[i, j] = bat_power[j]159 160 # 保存最佳粒子的位置161 pbest_all = np.zeros((sizepop, dim + 1))162 pbest = np.zeros((sizepop, 1))163 164 # 计算所有粒子的适应度函数值165 for i in range(0, sizepop):166 pbest[i] = fitness_caculate(pop_x[i, :])167 pbest_all[i, :] = pop_x[i, :]168 # 先默认第一个最好169 gbest_all = pop_x[0, :]170 gbest = pbest[0]171 172 # 全局最优更新,找出哪个粒子最优173 for i in range(1, sizepop):174 if pbest[i] < gbest:175 gbest_all[:] = pop_x[i, :]176 gbest = pbest[i]177 178 # 开始主循环179 uu = np.zeros(iteration)180 for t in range(iteration):181 for i in range(sizepop):182 w = w_max - (w_max - w_min) * t / iteration183 c1 = (0.5 - 2.5) * t / iteration + 2.5184 c2 = (2.5 - 0.5) * t / iteration + 0.5185 pop_v[i, :] = w * pop_v[i, :] + c1 * np.random.rand() * (186 pbest_all[i, :] - pop_x[i, :]) + c2 * np.random.rand() * (gbest_all - pop_x[i, :])187 for j in range(dim):188 if pop_v[i, j] > v2_max:189 pop_v[i, j] = v2_max190 elif pop_v[i, j] < v2_min:191 pop_v[i, j] = v2_min192 193 pop_x[i, :] += pop_v[i, :]194 195 for j in range(dim):196 if pop_x[i, j] > min(T_max - Load[j], 0.5 * pop_x[i, dim]):197 pop_x[i, j] = min(T_max - Load[j], 0.5 * pop_x[i, dim])198 pop_v[i, j] = -pop_v[i, j]199 elif (pop_x[i, j] + Load[j]) > T_max:200 pop_x[i, j] = T_max - Load[j]201 pop_v[i, j] = -pop_v[i, j]202 203 for j in range(dim):204 if pop_x[i, j] < max(-Load[j], -0.5 * pop_x[i, dim]):205 pop_x[i, j] = max(-Load[j], -0.5 * pop_x[i, dim])206 pop_v[i, j] = -pop_v[i, j]207 elif pop_x[i, j] < -(Load[j] + T_max):208 pop_x[i, j] = -(Load[j] + T_max)209 pop_v[i, j] = -pop_v[i, j]210 211 for j in range(dim):212 if pop_x[i, j] < 0 and Load[j] == 0:213 pop_x[i, j] = 0214 215 pop_x[i, dim] = ceil(pop_x[i, dim])216 217 y = fitness_caculate(pop_x[i, :])218 if y < pbest[i]:219 pbest_all[i, :] = pop_x[i, :]220 pbest[i] = y221 222 if pbest[i] < gbest:223 gbest_all[:] = pop_x[i, :]224 gbest = pbest[i]225 226 uu[t] = ceil(gbest)227 228 plt.figure(1)229 x1_ticks = np.arange(0, iteration + 1, 20)230 plt.xticks(x1_ticks, fontproperties='Times New Roman', size=16)231 plt.yticks(fontproperties='Times New Roman', size=16)232 plt.plot(uu)233 plt.xlabel('iteration', font={'family': 'Times New Roman', 'size': 16})234 plt.ylabel('Cost/RMB', font={'family': 'Times New Roman', 'size': 16})235 plt.rcParams['font.sans-serif'] = ['SimHei']236 plt.rcParams['axes.unicode_minus'] = False237 238 plt.figure(2)239 x2_ticks = np.arange(0, dim + 1, 80)240 y2_ticks = np.arange(-600, 200, 1000)241 242 plt.plot(gbest_all[0:dim], label='Energy storage power')243 plt.xticks(x2_ticks, ('0', '4', '8', '12', '16', '20', '24'), fontproperties='Times New Roman', size=16)244 plt.yticks(y2_ticks, fontproperties='Times New Roman', size=16)245 plt.plot(Load[0:dim], label='Load power')246 plt.xlabel('time/hour', font={'family': 'Times New Roman', 'size': 16})247 plt.ylabel('Power/kW', font={'family': 'Times New Roman', 'size': 16})248 # print('\n', '成本（元）', ceil(gbest), '\n', '储能容量（kWh）', gbest_all[dim], '\n')249 plt.legend()250 # plt.show()251 result = f'成本（元）: {ceil(gbest)};\n储能容量（kWh）: {gbest_all[dim]}.'252 pic_path = '.\\rf_Demo.png'253 plt.savefig(pic_path)254 return result, pic_path255 256 257 258 259demo = gr.Interface(260 fn = choose_file,261 inputs=["checkbox", gr.inputs.File()],262 outputs=["text", "image"],263 examples='',264 title="群体智能优化算法",265 description="请上传历史负荷电价数据或选择使用默认文件"266)267 268demo.launch()