因为在公司开发项目(新加坡),对于商用软件的使用较为敏感。所在项目需要使用到Cplex求解器,先利用python中可以直接使用的免费cplex库(限1000/1000)进行了简单的学习。
本文内容主要参考了此博客:Cplex +Python下的虚拟电厂智能调度实例(附代码)_cplex python 例子-CSDN博客
ps:因为AI的发展,现在的应用编程基本对代码的要求较低,初学可以利用一些AI软件进行辅助。
1. 借助conda软件(因为公司没买,这里我用的是miniconda)安装Python3.8,如果是vscode可以利用jupyter内核选择快速创建(点击右侧内核选择,然后点击中间选择其他内核,选择3.8版本)
2. 我这里conda库不能直接conda install,所以我是选定好环境后在notebook利用pip安装并重启内核。
pip install cplex docplex
- 以下是本次优化任务(初版来源于开头提到的参考博客)
虚拟电厂单用户调度(含储能)
- 光储荷经济性调度,针对单个用户、单个光伏、单个储能设备,短周期内的经济最优目标
变量
-
充电变量: X = ( x t ) ∈ [ 0 , α ] X = (x_t) \in [0, \alpha] X=(xt)∈[0,α]
-
放电变量: Y = ( y t ) ∈ [ − α , 0 ] Y = (y_t) \in [-\alpha, 0] Y=(yt)∈[−α,0]
-
储能设备能量状态: S = ( s t ) ∈ [ 0 , 1 ] S = (s_t) \in [0, 1] S=(st)∈[0,1]
以上三个变量为 continuous_var_list
连续变量列表
参数
-
负荷:
load
-
功率:
power
-
电价:
price
-
充电效率:
charge_efficiency
-
放电效率:
discharge_efficiency
-
储能单次充放电限制:
nominal_power
目标函数
最小化经济性最优目标:
minimize C o s t = ∑ t = 1 T ( l o a d ( t ) − p o w e r ( t ) + x ( t ) + y ( t ) ) × p r i c e ( t ) ÷ c o s t b a s e \text{minimize } Cost = \sum_{t=1}^T (load(t) - power(t) + x(t) + y(t)) \times price(t) \div costbase minimize Cost=t=1∑T(load(t)−power(t)+x(t)+y(t))×price(t)÷costbase
其中 costbase
为不使用储能的情况下的花费:
c
o
s
t
b
a
s
e
=
∑
t
=
1
T
(
l
o
a
d
(
t
)
−
p
o
w
e
r
(
t
)
)
×
p
r
i
c
e
(
t
)
costbase = \sum_{t=1}^T (load(t) - power(t)) \times price(t)
costbase=t=1∑T(load(t)−power(t))×price(t)
约束条件
-
储能能量状态:
-
储能初始状态为 0: s 0 = 0 s_0 = 0 s0=0
-
储能状态更新: s t = s t − 1 + c 1 ⋅ x t + y t c 2 s_t = s_{t-1} + c_1 \cdot x_t + \frac{y_t}{c_2} st=st−1+c1⋅xt+c2yt
-
from docplex.mp.model import Model
import random
#设置一个随机种子
random.seed(1234)
# 创建模型,在这里我们对连续10个状态进行动态优化,电价有4个时段
time = 4
time1 = 5
# 生成随机的用户负荷和用户功率
user_loads = [round(random.uniform(0, 1), 2) for _ in range(time1)]
user_powers = [round(random.uniform(0, 1), 2) for _ in range(time1)]
elec_price = [0.54, 0.22, 0.32, 0.24]
#设置储能充放电效率和单次充放限制
charge_eff = 0.91
discharge_eff = 0.95
nominal_power = 0.8
model = Model(name='Electricity')
x = model.continuous_var_list(time, lb=0,ub=nominal_power,name='x')
y = model.continuous_var_list(time, lb = -nominal_power,ub = 0, name='y')
soc = model.continuous_var_list(time1,lb=0,ub=nominal_power,name='soc')
#无储能下的花费
cost_base = sum(((user_loads[i] - user_powers[i]) * elec_price[i]) for i in range(time))
#总花费
total_cost = model.sum(((user_loads[i] - user_powers[i] + x[i] + y[i]) * elec_price[i] / cost_base) for i in range(time))
model.minimize(total_cost)
model.add_constraint(soc[0] == 0)
for i in range(time):
model.add_constraint(soc[i+1] == soc[i] + x[i] * charge_eff + y[i]/discharge_eff)
model.solve()
print('The optimized solution is:')
print("charging(+)/discharging(-) power for each time step:")
for i in range(time):
print("time{}:solution:".format(i),x[i].solution_value+y[i].solution_value)
for i in range(time1):
print('time{}:soc{}'.format(i,soc[i].solution_value))
print('Total cost:{}'.format(total_cost.solution_value))
4. 当我们把时间和状态扩展到24小时,每10分钟采样一次并指定策略,将电价设置为实时波动,可以将以上代码进行如下改进:
from docplex.mp.model import Model
import random
import matplotlib.pyplot as plt
# 设置一个随机种子
random.seed(1234)
# 创建模型,在这里我们对一天24小时进行动态优化,每10分钟采样一次
sampling_interval = 10 # 每10分钟采样一次
minutes_per_day = 24 * 60 # 一天共有 1440 分钟
time = minutes_per_day // sampling_interval # 采样点数量
# 生成随机的用户负荷和用户功率
user_loads = [round(random.uniform(0, 1), 2) for _ in range(time)]
user_powers = [round(random.uniform(0, 1), 2) for _ in range(time)]
# 生成随机电价,基于 0.5 生成 ±0.2 的波动
elec_price = [round(0.5 + random.uniform(-0.2, 0.2), 2) for _ in range(time)]
# 设置储能充放电效率和单次充放限制
charge_eff = 0.91
discharge_eff = 0.95
nominal_power = 0.8
# 创建模型
model = Model(name='Electricity')
# 创建变量
x = model.continuous_var_list(time, lb=0, ub=nominal_power, name='x')
y = model.continuous_var_list(time, lb=-nominal_power, ub=0, name='y')
soc = model.continuous_var_list(time, lb=0, ub=nominal_power, name='soc')
# 无储能下的花费
cost_base = sum(((user_loads[i] - user_powers[i]) * elec_price[i]) for i in range(time))
# 总花费
total_cost = model.sum(((user_loads[i] - user_powers[i] + x[i] + y[i]) * elec_price[i] / cost_base) for i in range(time))
# 设置目标函数
model.minimize(total_cost)
# 添加约束条件
model.add_constraint(soc[0] == 0)
for i in range(time - 1):
model.add_constraint(soc[i + 1] == soc[i] + x[i] * charge_eff + y[i] / discharge_eff)
# 求解模型
solution = model.solve()
# 检查是否找到最优解
if solution:
print('The optimized solution is:')
print("charging(+)/discharging(-) power for each time step:")
charging_discharge_values = []
soc_values = []
for i in range(time):
charge_discharge = x[i].solution_value + y[i].solution_value
charging_discharge_values.append(charge_discharge)
soc_values.append(soc[i].solution_value)
if i % (60 // sampling_interval) == 0: # 每小时输出一次状态
print(f"Hour {i // (60 // sampling_interval)}: solution: {charge_discharge}")
# 可视化结果
hours = [i * sampling_interval / 60 for i in range(time)]
plt.figure(figsize=(14, 10))
# 充放电功率可视化
plt.subplot(2, 1, 1)
plt.plot(hours, charging_discharge_values, label='Charging/Discharging Power', color='b')
plt.xlabel('Time (hours)')
plt.ylabel('Power (kW)')
plt.title('Charging (+) and Discharging (-) Power Over Time')
plt.legend()
# SOC 状态可视化
plt.subplot(2, 1, 2)
plt.plot(hours, soc_values, label='State of Charge (SOC)', color='g')
plt.xlabel('Time (hours)')
plt.ylabel('SOC')
plt.title('State of Charge Over Time')
plt.legend()
plt.tight_layout()
plt.show()
print(f'Total cost: {total_cost.solution_value}')
else:
print("No solution found.")
以上代码输出求解结果:
The optimized solution is:
charging(+)/discharging(-) power for each time step:
Hour 0: solution: 0.10839999999999994
Hour 1: solution: 0.2851999999999997
Hour 2: solution: 0.8
Hour 3: solution: 0.10840000000000005
Hour 4: solution: 0.8
Hour 5: solution: 0.8
Hour 6: solution: 0.8
Hour 7: solution: -0.7537304800462696
Hour 8: solution: 0.0
Hour 9: solution: -0.7537304800462696
Hour 10: solution: 0.0
Hour 11: solution: -0.6746096009253907
Hour 12: solution: 0.10839999999999994
Hour 13: solution: 0.17679999999999985
Hour 14: solution: 0.0
Hour 15: solution: 0.0
Hour 16: solution: 0.17679999999999985
Hour 17: solution: 0.10839999999999994
Hour 18: solution: -0.7537304800462696
Hour 19: solution: 0.0
Hour 20: solution: -0.6283400809716602
Hour 21: solution: 0.8
Hour 22: solution: 0.0
Hour 23: solution: 0.10839999999999994
我们可视化充放电状态和SOC状态:
- 此时我们能了解到单个系统内,一负载一发电光伏一储能并连接电网的简单模型优化,在此基础上长期优化还应当考虑以下约束:
- 增加复杂的约束条件:如电网容量、最小 SOC、电网平衡、冷却时间等。
- 改进目标函数:采用多目标优化,包括最小化总成本、最大化收益、平滑负荷等。
- 考虑不确定性和随机性:加入随机电价、负荷等因素,通过鲁棒优化提高模型的适用性。
- 引入实际电力市场的特性:如峰谷电价、惩罚机制等,确保模型能贴近实际运营情况。
- 设备寿命和维护:减少设备磨损,考虑储能充放电对设备寿命的影响。
- 比如我们这里可以考虑加入储能SOC下限(一般来说商用设备是在0.2~0.8间能够较好地延长使用寿命),并考虑到充放的冷却间隔(防止过热)
# 添加约束条件
min_soc = 0.2 * nominal_power # 最小 SOC 限制
model.add_constraint(soc[0] == 0.5 * nominal_power) # 初始 SOC 设置为 0.5
for i in range(time - 1):
# SOC 状态更新约束
model.add_constraint(soc[i + 1] == soc[i] + x[i] * charge_eff + y[i] / discharge_eff)
# 最小 SOC 约束
model.add_constraint(soc[i] >= min_soc)
# 放电后的冷却时间约束(通过引入额外的二进制变量来表示冷却状态)
for i in range(time - 2):
# 创建二进制变量来表示冷却状态
cooling_20 = model.binary_var(name=f'cooling_20_{i}')
cooling_10 = model.binary_var(name=f'cooling_10_{i}')
# 冷却 20 分钟约束,当放电量超过 50% 时触发冷却
model.add_indicator(cooling_20, y[i] <= -0.5 * nominal_power, 1)
model.add_constraint(x[i + 1] <= nominal_power * (1 - cooling_20)) # 冷却期间不可充电
# 冷却 10 分钟约束,当放电量在 30% - 50% 之间时触发冷却
model.add_indicator(cooling_10, y[i] <= -0.3 * nominal_power, 1)
model.add_constraint(x[i + 1] <= nominal_power * (1 - cooling_10)) # 冷却期间不可充电
可以看到这里的约束添加在随机数据中表现并不理想,可能原因很多,比如随机数据的波动让优化求解面临局限;约束的添加不合理等等。
以上内容均为本人个人实践和学习,不作为个人版权所有,欢迎大家指导讨论(可能存在很多潜在问题和研究方向,机器学习、联邦学习、能源优化等方向也欢迎邮件合作bo_chen@u.nus.edu(个人身份)