Python+Gurobi实现简单的微网优化运行小算例

之前一直用Matlab和工具箱的方式求解优化问题，今年暑期负责了暑期的Python课程，对Python有了更多的认识，Python确实是一个很简洁方便的语言。下面就以Python+Gurobi为例进行讲解。

准备工作

Python采用了Python3，Gurobi是9.1版本，Gurobi有academic license这个还蛮好的！开发环境为PyCharm。

如何在PyCharm环境使用Gurobi，可参考官网说明。
在程序中使用了numpy、scipy和matplotlib(一个绘图库)，需在Pycharm环境中通过pip install方式安装相应的库。

算例说明

在此假设一个并网型光储微电网，做一个24小时的优化运行问题，已知24小时负荷和光伏功率数据，在分时电价下，求解微电网的最优运行方案。
算例中的参数可根据自身情况做修改。

代码实现

import numpy as np
import scipy.sparse as sp
import gurobipy as gp
import matplotlib.pyplot as plt
from gurobipy import GRB

try:

    #定义输入数据，如负荷、光伏、电价信息等
    P_load = [700, 650, 500, 450, 600, 900, 1200, 1400, 1200, 1300, 1500, 1700, 1600, 1500, 1300, 1200, 1300, 1500, 1900, 1700, 1600, 1400, 1200, 1000] #24小时负荷数据
    P_pv = [0, 0, 0, 0, 0, 0, 200, 400, 800, 1200, 1500, 1800, 1950, 2000, 1800, 1500, 1000, 500, 100, 0, 0, 0, 0, 0] #24小时光伏数据

    Pricein = [0.43, 0.43, 0.43, 0.43, 0.43, 0.43, 0.43, 0.43, 0.69, 0.69, 0.69, 1.21, 1.21, 1.21, 1.21, 0.69, 0.69, 0.69, 1.21, 1.21, 1.21, 0.69, 0.69, 0.69] #购电电价
    Priceout = [0.27, 0.27, 0.27, 0.27, 0.27, 0.27, 0.27, 0.27, 0.5, 0.5, 0.5, 1.02, 1.02, 1.02, 1.02, 0.5, 0.5, 0.5, 1.02, 1.02, 1.02, 0.5, 0.5, 0.5] # 售电电价

    #定义系统相关参数，如电池参数，联络线功率参数等
    #电池参数
    Ess_max = 1500
    Ess_min = 300
    P_chmax = 500
    P_dismax = 500
    #联络线功率限制
    P_gridmax = 3000

    #定义优化时间
    T = 24

    #定义优化问题
    m = gp.Model("OptMG")

    #电池决策变量
    P_ch = m.addMVar(T, lb= 0.0, ub= P_chmax, vtype=GRB.CONTINUOUS, name="P_ch")
    P_dis = m.addMVar(T,lb=0.0, ub= P_dismax, vtype=GRB.CONTINUOUS, name="P_dis")
    F_ch = m.addMVar(T, vtype=GRB.BINARY, name="F_ch")
    F_dis = m.addMVar(T, vtype=GRB.BINARY, name="F_dis")
    Ess = m.addMVar(T+1, lb= Ess_min, ub=Ess_max, vtype=GRB.CONTINUOUS, name="P_dis")
    #电网交互功率决策变量
    P_gridin = m.addMVar(T, lb= 0.0, ub= P_gridmax, vtype=GRB.CONTINUOUS, name="P_gridin")
    P_gridout = m.addMVar(T, lb=0.0, ub=P_gridmax, vtype=GRB.CONTINUOUS, name="P_gridout")
    F_gridin = m.addMVar(T, vtype=GRB.BINARY, name="F_gridin")
    F_gridout = m.addMVar(T, vtype=GRB.BINARY, name="F_gridout")

    #目标函数
    obj = 0
    for t in range(T):
        obj += Pricein[t] * P_gridin[t] - Priceout[t] * P_gridout[t]

    m.setObjective(obj, GRB.MINIMIZE)

    #约束条件
    m.addConstr(Ess[0] == 300)
    m.addConstr(Ess[0] == Ess[T])
    for t in range(T):
        m.addConstr(P_load[t] + P_ch[t] + P_gridout[t] == P_pv[t] + P_dis[t] + P_gridin[t])
        m.addConstr(P_ch[t] <= P_chmax * F_ch[t])
        m.addConstr(P_dis[t] <= P_dismax * F_dis[t])
        m.addConstr(F_ch[t] + F_dis[t] <= 1)
        m.addConstr(P_gridin[t] <= P_gridmax * F_gridin[t])
        m.addConstr(P_gridout[t] <= P_gridmax * F_gridout[t])
        m.addConstr(F_gridin[t] + F_gridout[t] <= 1)

    for t in range(T-1):
        m.addConstr(Ess[t+1] == Ess[t] + P_ch[t] - P_dis[t])

    m.optimize()

    print('Obj: %g' % m.objVal)

    Pgrid = P_gridin.X - P_gridout.X
    P_ess = P_ch.X - P_dis.X
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    xlab = np.linspace(0,24,24)
    plt.plot(xlab,Pgrid,color='y',label=u'电网功率')
    plt.plot(xlab, P_ess, color='r',label=u'电池功率')
    plt.plot(xlab, P_load, color='b',label=u'负荷')
    plt.plot(xlab, P_pv, color='g',label=u'光伏')
    plt.show()

except gp.GurobiError as e:
    print('Error code ' + str(e.errno) + ': ' + str(e))

except AttributeError:
    print('Encountered an attribute error')