总结使用Pyomo解决优化问题的一般方式
首先当然要import pyomo.environ as pe
,以及定义m = pe.Concretemodel()
已知12个时刻的电价price_schedule
,以及12个时刻的充电量charge_schedule
求解目标是需要找到最好的售卖电量的方式 w t w_t wt使得总的利润 ∑ t = 1 T p t w t \sum_{t=1}^{\mathcal{T}}p_t w_t ∑t=1Tptwt最大,其中 p t p_t pt已知。
我在这里总结一下使用pyomo解答优化问题的关键要素。
1.定义参数的参数。这一部分包括总长度以及集合,规定了集合的总长度和形状。
m.nt = pe.Param(initialize=len(price_schedule), domain=pe.integers)
定义总长度之后,这只是一个数,print
以后可以得到12
,
用这个参数来定义时间集合T
,这是一个顺次生成的集合,其值为range(m.nt())
m.T = pe.Set(initialize=range(m.nt()))
有多少个总长度以及集合,上述过程就要重复多少次。当然,除外的情况为:这几个集合长度参数相同,因此只需要导入一个长度参数即可。但为了熟练多写几个长度参数也不是不可以。
2.定义参数。
在本问题中,包括 P t P_t Pt, Q t Q_t Qt, S 0 S_0 S0, w m a x w_{max} wmax四个参数。其中, P t P_t Pt, Q t Q_t Qt两个参数形状为m.T*1
,所以在pe.Param
中,应该加入对形状的说明,即m.T
即如下:
m.price = pe.Param(m.T, initialize=price_schedule)
m.charge = pe.Param(m.T, initialize=charge_schedule)
m.S0 = pe.Param(initialize=500.)
m.wmax = pe.Param(initialize=150.)
由于range(m.nt())
就是m.T
的值,所以可以range(m.nt())
用来指代m.T
。
3.定义变量。
这一段和往常相同,定义变量的内容使用pe.Var
即可。参数因为具有形状要求,应该加上range(m.nt())
或者m.T
m.w = pe.Var(m.T, domain=pe.NonNegativeReals)
m.s = pe.Var(m.T, domain=pe.NonNegativeReals)
4.定义目标
目标的表达式为
max w t , s t ∑ t ∈ T ( p t w t ) \begin{align} \max_{w_t, s_t}\quad& \sum_{t\in \mathcal{T}}(p_tw_t)\\ \end{align} wt,stmaxt∈T∑(ptwt)
目标为 max w t , s t ∑ t ∈ T ( p t w t ) \max_{w_t, s_t} \sum_{t\in \mathcal{T}}(p_tw_t) maxwt,st∑t∈T(ptwt),也就是 p t p_t pt和 w t w_t wt的对应项相乘。
原来的时候,目标使用m.o = pe.Objective(rule=lambda m:[optimization problem format], sense=maximize/minimize)
的形式定义。但是由于函数形式变长,这种定义形式变得繁琐,使用先定义优化函数,后定义目标的形式。
# 定义函数
def objective_func(model):return sum(m.price[t]*m.w[t] for t in m.T)
m.o = pe.Objective(rule=objective_func, sense=maximize)
5.定义约束
定义约束有
w t ≤ w m a x , ∀ t ∈ T S t ≤ S 0 , ∀ t ∈ T S t − S t − 1 = Q t − w t . ∀ t ∈ T , t ≥ 2 S t − S 0 = Q 1 − w 1 . \begin{align*} \quad & w_t\leq w_{max}, \forall t\in \mathcal{T}\\ & S_t\leq S^0 , \forall t\in \mathcal{T}\\ & S_{t}-S_{t-1}=Q_t-w_t. \forall t\in \mathcal{T},t\geq2\\ & S_{t}-S^0=Q_1-w_1. \end{align*} wt≤wmax,∀t∈TSt≤S0,∀t∈TSt−St−1=Qt−wt.∀t∈T,t≥2St−S0=Q1−w1.
在原来的时候,约束使用m.c = pe.ConstraintList()
形成约束列表,然后使用m.c.add()
的形式逐个写入约束。缺点是繁琐。
现在,先定义条件的函数def
,然后再使用model.[name of the def]=pe.Constriant()
导入约束条件。
- 可以用循环的形式理解这个代码,这里在函数调用过程中的
t
可以理解为循环单元,可以理解为类似于函数定义之类的东西,传入的是for t in T
里面的t
之类的东西。 - 在
pe.Constraint
中,等于对前面函数提到的t
进行迭代定义,从而得到一系列约束。再说一遍,m.T
即为range(m.nt())
,也就是从1
到m.nt()
的一个序列。
# 定义约束
def w_less_than_max(model, t):return m.w[t]<=m.wmax
m.w_less_than_max = pe.Constraint(m.T, rule=w_less_than_max)
类似地,定义第二个式子:
def Q_less_than_S0(model, t):return m.s[t]<=m.S0
m.Q_less_than_S0 = pe.Constraint(m.T, rule=Q_less_than_S0)
第三个式子有点复杂,分为2种情况,一种是t=1
的时候,一种则是其余时候
def Two_type_for_Storage(model,t):if t==0:return m.s[t]-m.s0==m.charge[t]*m.S0-m.w[t]else:return m.s[t]-m.s[t-1]==m.charge[t]*m.S0-m.w[t]
m.Two_type_for_Storage = pe.Constraint(m.T, rule=Two_type_for_Storage)
6.求解式子
和往常一样求解即可。
solver = pe.SolverFactory('gurobi')
results = solver.solve(model)
7.输出结果
# print(model.display())
print(f"time\tprice\tpower\tstorage")
for t in model.T:print(f"{t}\t{model.price[t]:.2f}\t{model.w[t]():>5.1f}\t{model.s[t]():>5.1f}")
也可以画图输出结果。这里就不赘述了。