数学建模

问题一

1、分析

假设火电机组1、2、3分别出力 $x_{1}MW$ 、 $x_{2}MW$ 、 $x_{3}MW$ ，则发电煤耗与其出力关系可以表示为 $F_{i} = a_{i}x_{i}^{2}+b_{i}x_{i}+c_{i}$ 。

则运行成本 = $\sum_{i=1}^{3}F_{i}*700/1000*1.5 = \sum_{i=1}^{3}F_{i}*1.05$

碳捕集成本 = $\sum_{i=1}^{3}CE_{i}*1000*x_{i}/1000*CC = \sum_{i=1}^{3}CE_{i}*x_{i}*CC$ ，其中 $CE_{i}$ 表示第 $i$ 个火点机组的碳排放量， $CC$ 表示单位碳捕集成本

2、建模

那么数学模型可以表示为：

min $\sum_{i=1}^{3}F_{i}*1.05+\sum_{i=1}^{3}CE_{i}*x_{i}*CC$

s.t. $\left\{\begin{matrix}C_{min}\le \sum_{i=1}^{3}x_{i}\le 900 \\180\le x_{1}\le 600\\90\le x_{2}\le 300\\45\le x_{3}\le 150\end{matrix}\right.$

3、求解

该问题是个二次规划问题，我们暂未采用商业求解器Gurobi、Cplex直接求解，也暂未使用启发式算法进行求解。

因此为了求解二次问题，有2种逼近思路

二次函数分段线性化：分段越多、结果越精准
决策变量设置步长前置计算二次函数值，决策变量从连续实数决策改为01决策：步长越短、结果越精准

目前我们采用了后者的求解思路，利用pyomo进行建模，调用开源求解器SCIP进行求解。

from pyomo.environ import *

power1, power2, power3 = list(), list(), list()
coal = 0

t = pd.read_excel("附件1.xlsx")
cmin_list = t["负荷功率(p.u.)"].to_list()

for r in cmin_list:

    # 创建一个模型实例
    model = ConcreteModel()

    def quadratic1(x):
        return 0.226*x**2+30.42*x+786.80
    def quadratic2(x):
        return 0.588*x**2+65.12*x+451.32
    def quadratic3(x):
        return 0.785*x**2+139.6*x+1049.50

    idx = list()
    for i in range(180, 601):
        idx.append((i, quadratic1(i)))
    model.x = Var(idx, domain=Binary)

    idx = list()
    for i in range(90, 301):
        idx.append((i, quadratic2(i)))
    model.y = Var(idx, domain=Binary)

    idx = list()
    for i in range(45, 151):
        idx.append((i, quadratic3(i)))
    model.z = Var(idx, domain=Binary)

    # 机组技术处理P
    p1 = sum([model.x[i, quadratic1(i)] * i for i in range(180, 601)])
    p2 = sum([model.y[i, quadratic2(i)] * i for i in range(90, 301)])
    p3 = sum([model.z[i, quadratic3(i)] * i for i in range(45, 151)])
    # 发电煤耗F
    f1 = sum([model.x[i, quadratic1(i)] * quadratic1(i) for i in range(180, 601)])
    f2 = sum([model.y[i, quadratic2(i)] * quadratic2(i) for i in range(90, 301)])
    f3 = sum([model.z[i, quadratic3(i)] * quadratic3(i) for i in range(45, 151)])
    # 定义目标
    cc=0
    model.obj = Objective(expr=((f1 + f2 + f3) * 1.05 + \
                          (0.72*p1 + 0.75*p2 + 0.79*p3) * cc), sense=minimize)

    # 定义约束
    model.constr1 = Constraint(expr=p1+p2+p3 <= 900)
    model.constr2 = Constraint(expr=900*r <= p1+p2+p3)

    model.constr3 = Constraint(expr=sum([model.x[i, quadratic1(i)] for i in range(180, 601)]) == 1)
    model.constr4 = Constraint(expr=sum([model.y[i, quadratic2(i)] for i in range(90, 301)]) == 1)
    model.constr5 = Constraint(expr=sum([model.z[i, quadratic3(i)] for i in range(45, 151)]) == 1)


    # 创建一个求解器实例
    solver = SolverFactory('scip')
    # model.write('test.lp')

    # 求解模型
    solver.solve(model)

    # # 打印结果
    # print("Solution:")
    # i, j, k = 0,0,0
    # for i in range(180, 601):
    #     if model.x[i, quadratic1(i)].value > 0:
    #         print(i, quadratic1(i))
    # for j in range(90, 301):
    #     if model.y[j, quadratic2(j)].value > 0:
    #         print(j, quadratic2(j))
    # for k in range(45, 151):
    #     if model.z[k, quadratic3(k)].value > 0:
    #         print(k, quadratic3(k))

    power1.append(sum([model.x[i, quadratic1(i)].value * i for i in range(180, 601)]))
    power2.append(sum([model.y[i, quadratic2(i)].value * i for i in range(90, 301)]))
    power3.append(sum([model.z[i, quadratic3(i)].value * i for i in range(45, 151)]))
    
    coal += sum([model.x[i, quadratic1(i)].value * quadratic1(i) for i in range(180, 601)]) + \
            sum([model.y[i, quadratic2(i)].value * quadratic2(i) for i in range(90, 301)]) + \
            sum([model.z[i, quadratic3(i)].value * quadratic3(i) for i in range(45, 151)])

4、结果

火电运行成本：coal*1.05/4/10000

碳捕集成本：(sum([0.72*power1[i] + 0.75*power2[i] + 0.79*power3[i] for i in range(96)]) * cc) / 4 /10000

总发电量（MWh）：(sum(power1)+sum(power2)+sum(power3)) / 4

碳捕集成本 (元/t)	火电运行成本（万元）	碳捕集成本（万元）	总发电成本（万元）	单位供电成本（万元/MWh）
0	244.886	0	244.886	0.016
60	244.894	68.331	313.225	0.020
80	244.899	91.102	336.001	0.022
100	244.907	113.869	358.776	0.023

图都基本一样，用碳捕集成本 = 0元/t，展示一个即可

问题二

1、分析

在问题一的模型结果基础上（注意：这里我们使用的是单位碳捕集成本为60元/t 的模型版本），如果风电装机300MW直接替代机组3

cloud = t["风电功率(p.u.)"].to_list()
print(len(cmin_list))
print(len([power1[i]+power2[i]+300*cloud[i] - 900*cmin_list[i]
for i in range(96) if power1[i]+power2[i]+300*cloud[i] > 900*cmin_list[i]]))
sum([power1[i]+power2[i]+300*cloud[i] - 900*cmin_list[i]
for i in range(96) if power1[i]+power2[i]+300*cloud[i] > 900*cmin_list[i]]) / 4

计算可以得到，此时系统96个时段的功率平衡中，会导致49个时段弃风，47个时段失电荷，其中弃风电量720.715MWh

为减少弃风又不失负荷，意思就是满足用户需求的前提下，最小化成本

2、建模

题干需要保证不失电荷，也就是必须满足用户需求，即不存在失电荷场景的成本

假设火电机组1、2，在t 时段分别出力 $x_{1, t}MW$ 、 $x_{2, t}MW$ ，风电 $x_{3}MW$ ，则发电煤耗与其出力关系可以表示为 $F_{i, t} = a_{i}x_{i, t}^{2}+b_{i}x_{i, t}+c_{i}$ 。

则全天火电成本 = $\sum_{i=1}^{2}\sum_{i=0}^{95}F_{i，t}*1.05+\sum_{i=1}^{2}\sum_{i=0}^{95}CE_{i}*x_{i,t}*CC$

全天风电成本 = $\sum_{t=0}^{96} 45*24*x_{3}*cloud[t]$

全天储能成本 = $\sum_{t=0}^{96}3000*1000*x_{3}*cloud[t]/10/365$

全天弃风成本 = $\sum_{i=0}^{95}(\sum_{i=1}^{2}x_{i,t}+x_{3}-900*cminlist[i])*300$

那么数学模型可以表示为：

min 上述成本相加

s.t. $\left\{\begin{matrix}\sum_{i=1}^{2}x_{i,t}+x_{3} \ge 900*cminlist[i],t\in range(96)\\180\le x_{1,t}\le 600,t\in range(96)\\90\le x_{2,t}\le 300,t\in range(96)\\0\le x_{3}\le 300\end{matrix}\right.$