水火电经济调度论文和程序.zip

#!/usr/bin/env python3 # -*- coding: utf-8 -*- # @Author: yudengwu(余登武） # @Date : 2021/4/15 #@email:1344732766@qq.com import numpy as np import pandas as pd from tqdm import tqdm # 进度条设置 import matplotlib.pyplot as plt from pylab import * from scipy import optimize as opt from scipy.optimize import minimize mpl.rcParams['font.sans-serif'] = ['SimHei'] mpl.rcParams['axes.unicode_minus'] = False """ 本例中考虑一个水库源头变化水火电系统的经济调度问题。该水火电系统包括两个火电站和两个水电站，且整个计划时间为24h， 其中，机组1和机组2为火电机组，机组3和机组4为水电机组。 具体机组数据见下： F1（PG1）=25+3.2*PG1+0.0025*PG1*PG1， 50MW<=PG1<=500MW F2(PG2)=30+3.4*PG2+0.0008*PG2*PG2 ，150MW<=PG2<=800MW 机组3（水利机组）的水流函数=0.1980+0.306*PG3+0.000216*PG3*PG3， 单位1e3 立方米/h 0<=PG3<=500MW 机组3（水力机组）的水库源头变化函数=0.90-0.0030*dk3+0.000001*dk3*dk3 机组4（水利机组）的水流函数=0.9360+0.612*PG4+0.000360*PG4*PG4 ，单位1e3 立方米/h 0<=PG4<=250MW 机组4（水力机组）的水库源头变化函数=0.95-0.0025*dk4+0.00002**dk4*dk4 两个水电站的水库表面积分别为1e7 和4e6 平方米，可用水资源分别为2.850e6 和2.450e6立方米，水库源头分别为30m和25m,整个计划期间水电站进水量为0 """ #水库源头可变情况下的输电损耗矩阵系数为 B=0.0001*np.array([[1.40,0.10,0.15,0.15], [0.15,0.60,0.10,0.13], [0.15,0.10,0.68,0.65], [0.15,0.13,0.65,0.70]]) # 读取数据 file = pd.read_csv('水火电经济调度数据.csv', encoding='gbk') PL=file['负荷需求/MW'] #非线性规划部分 def feixianxinguihua(PL): # 目标函数 def objective(X): # = B = 0.0001 * np.array([[1.40, 0.10, 0.15, 0.15], [0.15, 0.60, 0.10, 0.13], [0.15, 0.10, 0.68, 0.65], [0.15, 0.13, 0.65, 0.70]]) PG1 = X[0] PG2 = X[1] PG3 = X[2] PG4 = X[3] A = np.array([PG1, PG2, PG3, PG4]) A = A.reshape(1, 4) C = np.array([PG1, PG2, PG3, PG4]) C = C.reshape(4, 1) AB = np.dot(A, B) AC = np.dot(AB, C) ploss1t = AC + PG1 * B[0, 0] + PG2 * B[1, 0] + PG3 * B[2, 0] + PG4 * B[3, 0] + B[0, 0] ploss1t = np.squeeze(ploss1t) return (PG1 + PG2 + PG3 + PG4 - ploss1t - PL) # 约束条件 def constraint1(X): # 等式约束 B = 0.0001 * np.array([[1.40, 0.10, 0.15, 0.15], [0.15, 0.60, 0.10, 0.13], [0.15, 0.10, 0.68, 0.65], [0.15, 0.13, 0.65, 0.70]]) # 单位/MW PG1 = X[0] PG2 = X[1] PG3 = X[2] PG4 = X[3] A = np.array([PG1, PG2, PG3, PG4]) A = A.reshape(1, 4) C = np.array([PG1, PG2, PG3, PG4]) C = C.reshape(4, 1) AB = np.dot(A, B) AC = np.dot(AB, C) ploss1t = AC + PG1 * B[0, 0] + PG2 * B[1, 0] + PG3 * B[2, 0] + PG4 * B[3, 0] + B[0, 0] ploss1t = np.squeeze(ploss1t) return PG1 + PG2 + PG3 + PG4 - ploss1t - PL # 上下限 b1 = (50, 500) b2 = (150, 800) b3 = (0, 500) b4 = (0, 250) def run(): bnds = (b1, b2, b3, b4) con1 = {'type': 'eq', 'fun': constraint1} x0 = np.random.uniform(10, 800, 4) # 计算 solution = minimize(objective, x0, method='SLSQP', bounds=bnds, constraints=con1) if solution.success==True : x = solution.x return x[0], x[1], x[2], x[3] else:run() PG1,PG2,PG3,PG4=run() return PG1,PG2,PG3,PG4 ####################初始化参数##################### NP=50 #种群数量 L=5 # 对应5 个参数 4个机组出力 加功率损耗 Pc=0.5 #交叉率 Pm=0.1 #变异率 G=20 #最大遗传代数 t=24 #一天24小时 H_begin1=30 #水电站1 水库源头初始值 30m H_begin2=25 #水电站2 水库源头初始值 25m #H_end1=5 #水电站1 水库源头结束值 5m （这个值可以自己手动设置 #H_end2=5 #水电站2 水库源头结束值 5m （这个值可以自己手动设置 area1=1e7 #水库1垂直表面积 area2=4e6 #水库2垂直表面积 V1=2.850e6 #水库1可用水资源体积 V2=2.450e6 #水库2可用水资源体积 rh1=0 #水电站1进水速度 rh2=0 #水电站2进水速度 # 用一个 的矩阵表示种群，每行表示一个体 X = np.random.uniform(0, 500, size=(NP, 4, 24)) #机组1出力=X[,0,:],,机组2出力=X[：,1,:],机组3出力=X[;,2,:],机组4出力=X[；,3,:], ####################目标函数和惩罚项##################### #目标函数 #水火电经济调度 是通过所有在线运行的火力和水力机组进行分配，以使得在给定时间内满足功率平衡、发电机运行以及水力可用资源约束等情况下，火力发电机组的总发电成本最小。 #步骤2 目标函数 def calc_f(X): SUMCOST=[]#群体总成本 for i in range(NP): # 每一个个体 PG1=X[i,0,:]#机组1功率 火电机组 PG2=X[i,1,:]#机组2功率 火电机组 COST = [] # 存储个体费用 for j in range(t):#遍历每一个时刻 F=25 + 3.2 * PG1[j] + 0.0025 * PG1[j] * PG1[j]+30+3.4*PG2[j]+0.0008*PG2[j]*PG2[j] #当前时刻机组1和机组2成本 COST.append(F) SUMCOST.append(np.sum(COST)) return np.array(SUMCOST) #shape(50,) #步骤2 惩罚项 def calc_e(X): SUMCOST = [] # 存储群体惩罚项 for i in range(NP): # 每一个个体 PG1 = X[i, 0, :] # 机组1功率 火电机组 PG2 = X[i, 1, :] # 机组2功率 火电机组 PG3 = X[i, 2, :] #机组3功率 水电机组 PG4 = X[i, 3, :] #机组4功率 水电机组 COST = [] # 存储个体惩罚项 H3 = np.zeros((t+1, 1)) # 机组3的蓄水高度 H3[0,0]=30#初始高度 H4 = np.zeros((t + 1, 1)) # 机组4的蓄水高度 H4[0,0]=25#初始高度 # 2.蓄水约束 Q3T = [] # 存放机组3 每个时刻的水流速度 Q4T = [] # 存放机组4 每个时刻的水流速度 for j in range(t): # 遍历每一个时刻 q3t = (0.1980 + 0.306 * PG3[j] + 0.000216 * PG3[j] * PG3[j]) * ( 0.90 - 0.0030 * H3[j] + 0.000001 * H3[j] * H3[j]) # 机组3当前时刻水流速度。 q4t = (0.9360 + 0.612 * PG4[j] + 0.000360 * PG4[j] * PG4[j]) * ( 0.95 - 0.0025 * H4[j] + 0.00002 * H4[j] * H4[j]) # 机组4当前时刻水流速度 H3[j + 1] = H3[j] + 1e3*(rh1 - q3t) / area1 # 计算得到的下一时刻机组3水库源头高度 注意单位换算 要加1e3* H4[j + 1] = H4[j] + 1e3*(rh2 - q4t) / area2 # 计算得到的下一时刻机组3水库源头高度注意单位换算 要加1e3* Q3T.append(q3t) Q4T.append(q4t) if H3[j + 1] < 0: # 蓄水高度不能小于0 COST.append(np.abs(H3[j + 1])) if H4[j + 1] < 0: # 蓄水高度不能小于0 COST.append(np.abs(H4[j + 1])) # 在1个计划周期内 ，用水体积应小于等于可用水体积 if np.sum(Q3T) > V1/1000: # #print('np.sum(Q3T)',np.sum(Q3T)) COST.append(np.sum(Q3T) - V1/1000) if np.sum(Q4T) > V2/1000: #print('np.sum(Q4T)',np.sum(Q4T)) COST.append(np.sum(Q4