1 电力系统环境经济调度数学模型
2 改进粒子群算法解决
2.1 知识回顾
先回顾一下相关知识点：
2.2 案例1——IEEE6节点
我弄成了一个文档，对应起来：
（1）IEEE
#!/usr/bin/env python3 # -*- coding: utf-8 -*- # ~~~~~~~~~~~~~导入相关库~~~~~~~~~~· import numpy as np import pandas as pd import matplotlib.pyplot as plt import matplotlib as mpl mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体 mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题 import matplotlib matplotlib.use('TkAgg') from tqdm import tqdm # 进度条设置 # ~~~~~~~~~~~~~~~~~~~~~读取文件~~~~~~~~~~~~~~~~~~~~~~ file=pd.read_csv('IEEE6.csv')#机组参数 P_max = file['b'] #机组上限 P_min = file['c'] # #机组下限 # ~~~~~~~~~~~~~~~~~机组特性系数~~~~~~~~~~~~~~~~~~~~~~~~· ai = file['d'] bi=file['e'] ci=file['f'] # ~~~~~~~排放特性系数~~~~~~~~~~~~~~~~· ali=file['g'] beti=file['h'] gari=file['i'] # ~~~~~~~~~~~~~~~~~~~~~~~粒子群算法~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~· class Liqun: def __init__(self, Pload1, P_max,P_min,ai,bi,ci,ali,beti,gari): self.Pload1 = Pload1 # 负荷 self.P_max = P_max #机组上限 self.P_min = P_min #机组下限 self.ai=ai self.bi=bi self.ci=ci self.ali=ali self.beti=beti self.gari=gari # ~~~~~~~~PSO的参数~~~~~~~~~~~~~~~· self.w = 1 # 惯性因子，一般取1 self.c1 = 2 # 学习因子1，一般取2 self.c2 = 2 # 学习因子2，一般取2 self.m = 500 # 种群大小，即种群中小鸟的个数 self.iter_num = 100 # 算法最大迭代次数 self.max_vel = 0.5 # 限制粒子的最大速度为0.5 self.min_vel = -0.5 # 限制粒子的最小速度为-0.5 self.a1= 0.7 # 总燃料约束权重 self.a2= 0.3 # 污染排放量权重 # ~~~~~~~~~~~~~初始化群群体~~~~~~~~~~~~~~~· def InitializeX(self): """ :return: 群体[G1，G2，G3，G4，G5，G6] """ X = np.zeros((self.m, 6)) # 初始化群体，6代表 个机组出力 for n in range(self.m): # 遍历每一个粒子 X[n, 0] = np.random.uniform(10, 125, 1)[0] # G1 X[n, 1] = np.random.uniform(10, 150, 1)[0] # G2 X[n, 2] = np.random.uniform(35, 225, 1)[0] # G3 X[n, 3] = np.random.uniform(35, 210, 1)[0] # G4 X[n, 4] = np.random.uniform(130,325, 1)[0] # G5 X[n, 5] = np.random.uniform(125,315, 1)[0] # G6 return X # ~~~~~~~~~~~~~~~~~~~定义目标函数、和对应的惩罚项~~~~~~~~~~~~~~~~~~~· # ~~~~~1 目标函数：系统运行成本~~~~~ def function1(self, X1): """ 个体目标函数 :param X1: （个体[G1，G2，G3，G4，G5，G6] :return: 函数1目标函数值 """ SUMCOST = [] # 存储总的成本 for i in range (5): # 遍历每一台机组 cost = self.ci[i] * X1[i] * X1[i] + self.bi[i]* X1[i] +self.ai[i] # SUMCOST.append(cost) return np.sum(SUMCOST) #~~~~~~~~~~~~~~总污染排放量~~~~~~~~~~~~·· def function2(self, X1): """ 个体目标函数 :param X1: （个体[G1，G2，G3，G4，G5，G6] :return: 函数2目标函数值 """ emission=[] #储存总的污染排放量 for i in range(5): e=ali[i]+beti[i]*X1[i]+gari[i]*X1[i]*X1[i] emission.append(e) return np.sum(emission) # ~~~~平衡约束~~~~~~· def calc_e1(self, X1): """ 函数1 对应的个体惩罚项 :param X1: （个体[G1，G2，G3，G4，G5，G6] :return: """ cost=np.abs(X1[0]+X1[1]+X1[2]+X1[3]+X1[4]+X1[5]-self.Pload1) return cost # ~~~~~~~~~~~~~~~~粒子群速度更新公式~~~~~~~~~~~~~~~~~~~· def velocity_update(self, V, X, pbestX, gbestX): """ :param V: 群体速度 :param X: 群体位置[G1，G2，G3，G4，G5，G6] :param pbestX: 种群历史最优位置 :param gbestX: 全局最优位置 :return: """ r1 = np.random.random((self.m, 6)) r2 = np.random.random((self.m, 6)) V = self.w * V + self.c1 * r1 * (pbestX - X) + self.c2 * r2 * (gbestX - X) # 直接对照公式写就好了 # 防止越界处理 V[V > self.max_vel] = self.max_vel V[V < self.min_vel] = self.min_vel return V # ~~~~~~~粒子群位置更新公式~~~~~~~~~~~~~~~· def position_update(self, X, V): """ 根据公式更新粒子的位置 :param X: 粒子当前的位置矩阵，维度是 (self.m, 6) :param V: 粒子当前的速度举着，维度是 (self.m, 6) """ X = X + V # 更新位置 for n in range(self.m): # 遍历每一个粒子 if X[n, 0] < 10 or X[n, 0] > 125: # G1 X[n, 0] = np.random.uniform(10, 125, 1)[0] # G1 if X[n, 1] < 10 or X[n, 1] > 150: # G2 X[n, 1] = np.random.uniform(10, 150, 1)[0] # G2 if X[n, 2] < 35 or X[n, 2] > 225: # G3 X[n, 2] = np.random.uniform(35, 225, 1)[0] # G4 if X[n, 3] < 35 or X[n, 3] > 210: # G5 X[n, 3] = np.random.uniform(35, 210, 1)[0] # G8 if X[n, 4] < 130 or X[n, 4] > 325: # G11 X[n, 4] = np.random.uniform(130, 325, 1)[0] # G1 X[n, 5] = self.Pload1 - X[n, 0] - X[n, 1] - X[n, 2] - X[n, 3] - X[ n, 4] # 负荷平衡得到 G13 出力 while (X[n, 5] < 125) or (X[n, 5] > 315): # 如果G13 上下限违规 X[n, 0] = np.random.uniform(10, 125, 1)[0] # G1 X[n, 1] = np.random.uniform(10, 150, 1)[0] # G2 X[n, 2] = np.random.uniform(35, 225, 1)[0] # G3 X[n, 3] = np.random.uniform(35, 210, 1)[0] # G4 X[n, 4] = np.random.uniform(130, 325, 1)[0] # G5 X[n, 5] = self.Pload1- X[n, 0] - X[n, 1] - X[n, 2] - X[n, 3] - X[ n, 4] # 负荷平衡得到 G13 出力 return X
（2）结果
100%|██████████| 100/100 [00:03<00:00, 26.09it/s] 迭代最优变量是： [ 41.42608737 36.4634904 182.18059225 134.44146782 290.55563308 314.93272908] 迭代最优目标函数1(总燃料成本）：35260.71417 迭代最优目标函数2(总污染排放量）：605.34109 迭代最优的惩罚项： 0.0
2.3 案例2——IEEE10
接下来是10机组有阀点效应的，有网损的。最后再来2.4节40机组的，难度慢慢加大、第3节额外加了点餐：
多目标遗传算法NSGAII求解微电网调度(风、光、电动车、柴油机、主网)【python源码实现】
（1）Python代码
# ~~~~~~~~~~~~~~~~~~~~~~~粒子群算法~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~· class Liqun: def __init__(self, Pload1, P_max,P_min,ai,bi,ci,di,ei,ali,beti,gari,eti,dali): self.Pload1 = Pload1 # 负荷 self.P_max = P_max #机组上限 self.P_min = P_min #机组下限 self.ai=ai self.bi=bi self.ci=ci self.di=di self.ei=ei self.ali=ali self.beti=beti self.gari=gari self.eti=eti self.dali=dali # ~~~~~~~~PSO的参数~~~~~~~~~~~~~~~· self.w = 1 # 惯性因子，一般取1 self.c1 = 2 # 学习因子1，一般取2 self.c2 = 2 # 学习因子2，一般取2 self.m = 500 # 种群大小，即种群中小鸟的个数 self.iter_num = 100 # 算法最大迭代次数 self.max_vel = 0.5 # 限制粒子的最大速度为0.5 self.min_vel = -0.5 # 限制粒子的最小速度为-0.5 self.a1= 0.7 # 总燃料约束权重 self.a2= 0.3 # 污染排放量权重 # ~~~~~~~~~~~~~初始化群群体~~~~~~~~~~~~~~~· def InitializeX(self): """ :return: 群体[G1，G2，G3，G4，G5，G6,G7,G8,G9,G10] """ X = np.zeros((self.m, 10)) # 初始化群体，6代表 个机组出力 for n in range(self.m): # 遍历每一个粒子 X[n, 0] = np.random.uniform(10, 55, 1)[0] # G1 X[n, 1] = np.random.uniform(20, 80, 1)[0] # G2 X[n, 2] = np.random.uniform(47, 120, 1)[0] # G3 X[n, 3] = np.random.uniform(20, 130, 1)[0] # G4 X[n, 4] = np.random.uniform(50,160, 1)[0] # G5 X[n, 5] = np.random.uniform(70,240, 1)[0] # G6 X[n, 6] = np.random.uniform(60, 300, 1)[0] # G7 X[n, 7] = np.random.uniform(70, 340, 1)[0] # G8 X[n, 8] = np.random.uniform(130, 470, 1)[0] # G9 X[n, 9] = np.random.uniform(150, 470, 1)[0] # G10 return X # ~~~~~~~~~~~~~~~~~~~定义目标函数、和对应的惩罚项~~~~~~~~~~~~~~~~~~~· # ~~~~~1 目标函数：系统运行成本~~~~~ def function1(self, X1): """ 个体目标函数 :param X1: （个体[G1，G2，G3，G4，G5，G6,G7,G8,G9,G10] :return: 函数1目标函数值 """ SUMCOST = [] # 存储总的成本 for i in range (9): # 遍历每一台机组,Python是从0开始索引，0-9：10台机组 cost = self.ci[i] * X1[i] * X1[i] + self.bi[i]* X1[i] +self.ai[i]+\ np.abs(self.di[i]*np.sin(ei[i]*(self.P_min[i]-X1[i]))) #考虑阀点效应 SUMCOST.append(cost) return np.sum(SUMCOST) #~~~~~~~~~~~~~~总污染排放量~~~~~~~~~~~~·· def function2(self, X1): """ 个体目标函数 :param X1: （个体[G1，G2，G3，G4，G5，G6,G7,G8,G9,G10] :return: 函数2目标函数值 """ emission=[] #储存总的污染排放量 for i in range(9): e=self.ali[i]+self.beti[i]*X1[i]+self.gari[i]*X1[i]*X1[i]+self.eti[i]*np.exp(self.dali[i]*X1[i]) emission.append(e) return np.sum(emission) # ~~~~平衡约束~~~~~~· def calc_e1(self, X1): """ 函数1 对应的个体惩罚项 :param X1: （个体[G1，G2，G3，G4，G5，G6,G7,G8,G9,G10] :return: """ Ploss=0 for i in range(9): for j in range(9): Ploss+=X1[i]*B[i][j]*X1[j] # cost=np.abs(X1[0]+X1[1]+X1[2]+X1[3]+X1[4]+X1[5]+X1[6]+X1[7]+X1[8]+X1[9]-np.sum(Ploss)-self.Pload1) # return cost return max(0,abs(sum(X1)-Ploss-Pload1-0.1))
（2）结果
100%|██████████| 100/100 [00:15<00:00, 6.64it/s] 迭代最优变量是： [ 37.64814287 22.49213541 47.1752439 28.80811605 71.07470515 163.16909614 103.42579949 162.7489263 158.39058084 311.7171866 ] 迭代最优目标函数1(总燃料成本）：46726.82336 迭代最优目标函数2(总污染排放量）：4708365.86991 Process finished with exit code 0
2.4 案例3——IEEE40
（1）代码
#!/usr/bin/env python3 # -*- coding: utf-8 -*- # ~~~~~~~~~~~~~导入相关库~~~~~~~~~~· import numpy as np import pandas as pd import matplotlib.pyplot as plt import matplotlib as mpl mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体 mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题 import matplotlib matplotlib.use('TkAgg') from tqdm import tqdm # 进度条设置 # ~~~~~~~~~~~~~~~~~~~~~读取文件~~~~~~~~~~~~~~~~~~~~~~ file=pd.read_csv('IEEE40.csv')#机组参数 P_max = file['b'] #机组上限 P_min = file['c'] # #机组下限 # ~~~~~~~~~~~~~~~~~机组特性系数~~~~~~~~~~~~~~~~~~~~~~~~· ai = file['d'] bi=file['e'] ci=file['f'] di=file['g'] ei=file['h'] # ~~~~~~~排放特性系数~~~~~~~~~~~~~~~~· ali=file['g'] beti=file['h'] gari=file['i'] eti=file['l'] dali=file['m'] # ~~~~~~~~~~~~~~~~~~~~~~~粒子群算法~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~· class Liqun: def __init__(self, Pload1, P_max,P_min,ai,bi,ci,di,ei,ali,beti,gari,eti,dali): self.Pload1 = Pload1 # 负荷 self.P_max = P_max #机组上限 self.P_min = P_min #机组下限 self.ai = ai self.bi = bi self.ci = ci self.di = di self.ei = ei self.ali = ali self.beti = beti self.gari = gari self.eti = eti self.dali = dali # ~~~~~~~~PSO的参数~~~~~~~~~~~~~~~· self.w = 1 # 惯性因子，一般取1 self.c1 = 2 # 学习因子1，一般取2 self.c2 = 2 # 学习因子2，一般取2 self.m = 500 # 种群大小，即种群中小鸟的个数 self.iter_num = 100 # 算法最大迭代次数,机组数多，可以迭代次数加大 self.max_vel = 0.5 # 限制粒子的最大速度为0.5 self.min_vel = -0.5 # 限制粒子的最小速度为-0.5 self.a1= 0.4 # 总燃料约束权重 self.a2= 0.6 # 污染排放量权重
(2)结果
100%|██████████| 100/100 [00:30<00:00, 3.33it/s] 迭代最优变量是： [ 30.34869396 51.10965795 98.71908161 87.97301056 70.12927254 112.97742184 109.20125594 145.88263194 140.71615841 198.74327982 315.80318374 301.38126233 228.79039464 265.50721731 473.09954507 196.12949192 428.10579218 455.2538421 484.26009872 287.85042829 500.34306436 533.88884883 526.00442712 449.67581489 359.71492346 406.90418265 102.02746082 115.54292506 95.39280101 76.10125624 98.90489259 137.27482259 148.20881091 165.974008 178.99502653 165.49481492 70.68453874 48.12768985 69.87007254 1.87005299] 迭代最优目标函数1(总燃料成本）：9975.28353 迭代最优目标函数2(总污染排放量）：10314225.31643 Process finished with exit code 0
3 Python代码
# ~~~~~~~~~~~~~~~~~~~~~~~粒子群算法~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~· class Liqun: def __init__(self, Pload1, P_max,P_min,ai,bi,ci,di,ei,ali,beti,gari,eti,dali): self.Pload1 = Pload1 # 负荷 self.P_max = P_max #机组上限 self.P_min = P_min #机组下限 self.ai=ai self.bi=bi self.ci=ci self.di=di self.ei=ei self.ali=ali self.beti=beti self.gari=gari self.eti=eti self.dali=dali # ~~~~~~~~PSO的参数~~~~~~~~~~~~~~~· self.w = 1 # 惯性因子，一般取1 self.c1 = 2 # 学习因子1，一般取2 self.c2 = 2 # 学习因子2，一般取2 self.m = 500 # 种群大小，即种群中小鸟的个数 self.iter_num = 100 # 算法最大迭代次数 self.max_vel = 0.5 # 限制粒子的最大速度为0.5 self.min_vel = -0.5 # 限制粒子的最小速度为-0.5 self.a1= 0.7 # 总燃料约束权重 self.a2= 0.3 # 污染排放量权重 # ~~~~~~~~~~~~~初始化群群体~~~~~~~~~~~~~~~· def InitializeX(self): """ :return: 群体[G1，G2，G3，G4，G5，G6,G7,G8,G9,G10] """ X = np.zeros((self.m, 10)) # 初始化群体，6代表 个机组出力 for n in range(self.m): # 遍历每一个粒子 X[n, 0] = np.random.uniform(10, 55, 1)[0] # G1 X[n, 1] = np.random.uniform(20, 80, 1)[0] # G2 X[n, 2] = np.random.uniform(47, 120, 1)[0] # G3 X[n, 3] = np.random.uniform(20, 130, 1)[0] # G4 X[n, 4] = np.random.uniform(50,160, 1)[0] # G5 X[n, 5] = np.random.uniform(70,240, 1)[0] # G6 X[n, 6] = np.random.uniform(60, 300, 1)[0] # G7 X[n, 7] = np.random.uniform(70, 340, 1)[0] # G8 X[n, 8] = np.random.uniform(130, 470, 1)[0] # G9 X[n, 9] = np.random.uniform(150, 470, 1)[0] # G10 return X # ~~~~~~~~~~~~~~~~~~~定义目标函数、和对应的惩罚项~~~~~~~~~~~~~~~~~~~· # ~~~~~1 目标函数：系统运行成本~~~~~ def function1(self, X1): """ 个体目标函数 :param X1: （个体[G1，G2，G3，G4，G5，G6,G7,G8,G9,G10] :return: 函数1目标函数值 """ SUMCOST = [] # 存储总的成本 for i in range (9): # 遍历每一台机组,Python是从0开始索引，0-9：10台机组 cost = self.ci[i] * X1[i] * X1[i] + self.bi[i]* X1[i] +self.ai[i]+\ np.abs(self.di[i]*np.sin(ei[i]*(self.P_min[i]-X1[i]))) #考虑阀点效应 SUMCOST.append(cost) return np.sum(SUMCOST) #~~~~~~~~~~~~~~总污染排放量~~~~~~~~~~~~·· def function2(self, X1): """ 个体目标函数 :param X1: （个体[G1，G2，G3，G4，G5，G6,G7,G8,G9,G10] :return: 函数2目标函数值 """ emission=[] #储存总的污染排放量 for i in range(9): e=self.ali[i]+self.beti[i]*X1[i]+self.gari[i]*X1[i]*X1[i]+self.eti[i]*np.exp(self.dali[i]*X1[i]) emission.append(e) return np.sum(emission) # ~~~~平衡约束~~~~~~· def calc_e1(self, X1): """ 函数1 对应的个体惩罚项 :param X1: （个体[G1，G2，G3，G4，G5，G6,G7,G8,G9,G10] :return: """ Ploss=0 for i in range(9): for j in range(9): Ploss+=X1[i]*B[i][j]*X1[j] # cost=np.abs(X1[0]+X1[1]+X1[2]+X1[3]+X1[4]+X1[5]+X1[6]+X1[7]+X1[8]+X1[9]-np.sum(Ploss)-self.Pload1) # return cost return max(0,abs(sum(X1)-Ploss-Pload1-0.1))
