問題一
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，展示即可
