总结使用Pyomo解决优化问题的一般方式

总结使用Pyomo解决优化问题的一般方式

首先当然要import pyomo.environ as pe，以及定义m = pe.Concretemodel()

已知12个时刻的电价price_schedule,以及12个时刻的充电量charge_schedule

求解目标是需要找到最好的售卖电量的方式$w_t$使得总的利润${\sum }_{t=1}^{\mathcal{T}}{p}_t{w}_t$最大，其中$p_t$已知。

我在这里总结一下使用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$，$Q_t$，$S_0$，$w_{max}$四个参数。其中，$P_t$，$Q_t$两个参数形状为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.定义目标

目标的表达式为
$$\begin{array}{rl}\underset{w_t,s_t}{\max }& \sum _{t\in \mathcal{T}}(p_t{w}_t)\end{array}$$
目标为${\max }_{w_t,s_t}{\sum }_{t\in \mathcal{T}}(p_t{w}_t)$，也就是$p_t$和$w_t$的对应项相乘。

原来的时候，目标使用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.定义约束

定义约束有
$$\begin{array}{rl}& w_t\le w_{max},\forall t\in \mathcal{T}\\ & S_t\le S^0,\forall t\in \mathcal{T}\\ & S_t-S_{t-1}=Q_t-w_t.\forall t\in \mathcal{T},t\ge 2\\ & S_t-S^0=Q_1-w_1.\end{array}$$
在原来的时候，约束使用m.c = pe.ConstraintList()形成约束列表，然后使用m.c.add()的形式逐个写入约束。缺点是繁琐。

现在，先定义条件的函数def，然后再使用model.[name of the def]=pe.Constriant()导入约束条件。

1. 可以用循环的形式理解这个代码，这里在函数调用过程中的t可以理解为循环单元，可以理解为类似于函数定义之类的东西，传入的是for t in T里面的t之类的东西。
2. 在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}")

也可以画图输出结果。这里就不赘述了。