找到Pyomo最优解，但解数为0（水电模拟）

Question

我希望模拟挪威的水力发电，考虑到所有的梯级依赖性（每小时一次的解决方案）。首先，我在Pyomo设置了一个优化问题，针对一天中的三个工厂和三个水库。目前，我想最大限度地提高所有时间和所有发电厂/水库的利润，计算公式如下：电价[$/MWh]*水轮机出力[MW]*3600[s=1h]+容积[m3]*能源当量[MWh/m3]*电价[$/MWh]

我的限制是：

小时负荷=发电量必须等于每小时的负荷

水平衡=在每个时间步正确更新存储水平

开始和结束蓄水位=模拟期开始和结束时最大蓄水量的60%（这处于水平衡限制条件下）

问题：

1）为什么print语句为cascade_inflow和excleded_flow打印0.0：因为这是在创建模型时打印的，而不是在求解模型时打印的？在

2A）终止条件是最优的，我对目标函数有一个值，但解的数目是0：我认为问题在于约束水平衡，我将把结果再发布6小时。我是否需要设置水平衡的上限和下限限制？ 编辑：如果我从命令行使用pyomo solve--solver=glpk运行脚本脚本.py输入.dat我得到了一个解决方案。。？！在

2B）水约束功能不正常。如果我看一下结果，那么从时间步骤1到2的数量跳跃是不可行的。怎么了？我添加级联流的方式是不正确的，还是变量m.volume只是做它想要的？在

3）创建网络流问题是否是一个更好的主意？在Pyomo Gallery中有这类问题的示例代码。但我还不确定如何将节点建模为水库。（我很可能会在尝试实现一个脚本后立即为此发布一个新的帖子）。在

（这是我的第一篇博文：我做错了什么或者下次应该做得更好吗？）在

代码（读取忽略的参数）

from pyomo.environ import *
from pyomo.opt import SolverFactory

opt = SolverFactory("glpk")

# Initiate model instance
m = AbstractModel()

# Define the variable power for each time step
m.turbine = Var(m.stage, m.res, initialize=0, bounds=turbine_bounds, within=NonNegativeReals)
# Define the variable volume for each time step
m.volume = Var(m.stage, m.res, initialize=volume_Init,  bounds=volume_bounds, within=NonNegativeReals)
# Define the variable spill flow for each time step
m.spilledFlow = Var(m.stage, m.res, initialize=0, bounds=spill_bounds, within=NonNegativeReals)


# Constrain total power generated over day
def hourly_load_rule(model, t):
    return model.P_load*model.hourly_demand[t] <= sum(model.turbine[t, res] for res in model.res) <= model.P_load*model.hourly_demand[t]
m.hourly_load = Constraint(m.stage, rule=hourly_load_rule)


# Water balance equation
def water_balance_rule(model, t, r):
    if t == model.T:  # final volume is same as initial volume at 60%
        return model.volume[t, r] == model.max_Vol[r]*0.6
    elif t > 1:
        cascade_inflow = 0
        for stor in model.res:
            # connectMat is a matrix that has a 1 where there is a connection and 0 where there is not
            cascade_inflow = cascade_inflow + model.connectMat[stor, r]*(model.turbine[t, stor]/model.slope[stor]+model.spilledFlow[t, stor])
            if model.connectMat[stor, r] == 1:
                print(stor, r, t, value(cascade_inflow))  # this always prints out 0.0 for cascade_inflow
         discharged_flow = model.turbine[t, r]/model.slope[r]  # model.turbine is in MW: divide by slope to get discharge flow [m3/s]
         print(value(discharged_flow))  # this always prints out 0.0
         return model.volume[t, r] == model.volume[t-1, r]+(cascade_inflow+model.inflow[t, r]-model.spilledFlow[t, r]-discharged_flow)*model.secPerTimeStep
    else:
        # start volume constrained to 60% of max volume
        return model.volume[t, r] == model.max_Vol[r]*0.6
m.water_balance = Constraint(m.stage, m.res, rule=water_balance_rule)


# Revenue = Objective function (sum over all hours and all plants/reservoirs)
def revenue_rule(model):
    return sum(sum(model.el_price[i]*model.turbine[i, j]*model.secPerTimeStep+model.volume[i, j]*model.energy_stored[j]*model.el_price[i] for i in model.stage) for j in model.res)
m.obj = Objective(rule=revenue_rule, sense=maximize)

instance = m.create("three_input_stack.dat")

results = opt.solve(instance)

instance.display()

results.write()

结果

^{pr2}$

输入文件

param secPerTimeStep:=3600;
param T:=6;
param numReservoirs:=3;
param connectMat:=
1   1   0
1   2   1
1   3   0
2   1   0
2   2   0
2   3   1
3   1   0
3   2   0
3   3   0;
param el_price:=
1   242.16
2   242.09
3   239.3
4   231.52
5   224.25
6   219.77;
param inflow:=
1   1   5
2   1   5
3   1   5
4   1   5
5   1   5
6   1   5
1   2   5
2   2   5
3   2   5
4   2   5
5   2   5
6   2   5
1   3   5
2   3   5
3   3   5
4   3   5
5   3   5
6   3   5;
param min_Vol:=
1   0.0
2   0.0
3   0.0;
param max_Vol:=
1   66000000.0
2   16700000.0
3   2100000.0;
param min_Turb_gen:=
1   0
2   0
3   0;
param max_Turb_gen:=
1   3.31
2   5.9
3   9.0;
param min_spill:=
1   0.0
2   0.0
3   0.0;
param max_spill:=
1   10000.0
2   10000.0
3   10000.0;
param min_discharge:=
1   0.0
2   0.0
3   0.0;
param max_discharge:=
1   20.0
2   15.0
3   8.0;
param slope:=
1   0.1655
2   0.3933
3   1.125;
param energy_stored:=
1   0.000046
2   0.000109
3   0.00031;
param hourly_demand:=
1   0.0351
2   0.0335
3   0.0328
4   0.0325
5   0.0331
6   0.0357;
param P_load:=200;

Answer 1

1）总是打印零，因为m.tubine被初始化为0。在模型构造期间（即执行print语句时），表达式的值为0。在

2A）调用解算器后，结果会自动加载回模型中，因此命令instance.display()正在打印结果。results对象仅在您想延迟加载解决方案时使用（这是在调用解算器时可以设置的选项）。在

2B）不确定。您可以始终运行命令instance.water_balance.pprint()，以确认约束表达式是您所期望的。在