对多个位置的PYOMO小时优化：为什么所有输出变量等于零？

Question

简短版本

我正在开发一种优化算法，以最大程度地减少成本函数。由于我想将此功能应用于多个独立的位置，因此将最小化应用于所有位置的成本函数的 sum （因此，我不必为每个位置应用一个目标函数）。此外，成本计算是年度，但始于小时值。因此，对于每个位置，必须将所有小时的值总结。
该模型在下面详细介绍。

我的问题是：当我运行程序时，所有变量仅获得零；自然，这不是我想要的解决方案，但它甚至不应该是一个可能的解决方案 - 因为目标功能是商。
另外，求解器表示，当我在一个位置运行时，有35043个变量 + 35041的约束，而175215变量 + 175205我尝试五个位置时的约束 - 看起来太多了；对我来说真的很奇怪。

据我所知，我已经尽力编写代码的替代版本（更好，更轻松的）版本：消除冗余支持变量，将平等约束转换为不平等约束，并将目标函数从非线性转换为二次。似乎没有任何改变结果。

长期

摘要

我旨在建立一种优化算法，以找到在给定的一组合适位置中生产1 kg绿色氢的最低成本。

每个位置中可能存在两个发电厂 - 太阳能光伏公园和一个风电场 - 以及一个电解室系统。发电厂的安装容量和电解仪的标称输入功率是决策变量。也就是说，我想知道每个位置的这些变量的值是什么。为此，我必须对小时值进行一系列操作，以将电力转换为氢，然后评估一年内的总产量。我认为这是问题所在的地方 - 连接位置和小时的集合。

下面描述的所有参数是每个位置（和每个小时）的标量相同的标量，除了容量因子，每小时会变化并索引到每个位置。

由于这是一个确定性优化，其中每个位置独立于相邻的优化，因此可以进行所有合适位置的联合最小化（总和的最小化等于最小化的总和）。

问题符号

可以使用以下模型以数学形式以数学形式提出该问题：

模型初始化

from __future__ import division
from pyomo.opt import SolverFactory
from pyomo.core import AbstractModel
from pyomo.dataportal.DataPortal import DataPortal
from pyomo.environ import *
import pandas as pd
import os
os.environ["octeract_options"] = "num_cores=1"

model = AbstractModel()
data = DataPortal()

集

从预先准备的。csv 文件，可用
l ：合格的位置点
h ：一年

model.L = Set()
data.load(filename='modelL-5.csv', set=model.L)

model.H = Set()
data.load(filename='modelH.csv', set=model.H)

参数

的小时 k_solar ：太阳能公园总成本
k_wind ：风电场总成本
k_h2 ：电解器总成本

pmax_solar ：太阳能公园最大安装容量
pmax_wind ：风电场最大安装容量

pden_solar ：太阳能公园电源密度
pden_wind ：风电源电源密度
pden_h2 ：电解器功率密度

alpha ：每个位置的可用面积

effi ：electrolyser效率

sy ：生产扩散因子

K_solar = 728.5
K_wind = 1559.5
K_h2 = 1063.9

Pmax_solar = 50300
Pmax_wind = 475960

Pden_solar = 0.08268
Pden_wind = 0.78231
Pden_h2 = 0.01172

alpha = 608400

effi = 0.0159

Sy = 11.9

< em>容量因子索引到每个位置和每个小时；从预先准备的。csv files /s！alqcyawzssicg6aseth7u0edilm78a？ CFSolar ：太阳能因子
cfwind ：风能因素

model.CFsolar = Param(model.L, model.H)
data.load(filename='CF_Solar-5.csv', param=model.CFsolar, format='array')

model.CFwind = Param(model.L, model.H)
data.load(filename='CF_Wind-5.csv', param=model.CFwind, format='array')

决策变量

（我真正想要的）
仅取决于位置的变量：
psolar ：太阳能公园的安装能力
pwind ：风电场的安装能力
ph2 ：电解器额定功率

model.Psolar = Var(model.L, initialize=5000, domain=NonNegativeReals, bounds=(0, Pmax_solar))
model.Pwind = Var(model.L, initialize=5000, domain=NonNegativeReals, bounds=(0, Pmax_wind))
model.Ph2 = Var(model.L, initialize=5000, domain=NonNegativeReals, bounds=(0, Pmax_solar+Pmax_wind))

支撑变量

（我不想要的功率，但需要获取决策变量）
取决于一年的位置和小时
的变量 e_res ：产生的电力
e_h2 ：电力
的电力 e_curt ：电力减少
y_h2 ：由电解器

model.E_res = Var(model.L, model.H, domain=NonNegativeReals, bounds=(0, Pmax_solar+Pmax_wind))
model.E_h2 = Var(model.L, model.H, domain=NonNegativeReals, bounds=(0, Pmax_solar+Pmax_wind))
model.E_curt = Var(model.L, model.H, domain=NonNegativeReals, bounds=(0, Pmax_solar+Pmax_wind))
model.Y_h2 = Var(model.L, model.H, domain=NonNegativeReals)

数学公式

约束方程

产生的氢，每个位置的电力小时产生等于该位置上的太阳能/风能安装的电源相应的容量因子（取决于位置的位置因素和小时）

def RESenergy_production_rule(model, l, h):
    return model.E_res[l, h] == model.Psolar[l] * model.CFsolar[l, h] + model.Pwind[l] * model.CFwind[l, h]
model.RESenergy_production = Constraint(model.L, model.H, rule=RESenergy_production_rule)

对于每个位置，在每个小时内，该产生的电能都可以直接到电解器或未使用（限制）

def RESenergy_use_rule(model, l, h):
    return model.E_res[l, h] == model.E_h2[l, h] + model.E_curt[l, h]
model.RESenergy_use = Constraint(model.L, model.H, rule=RESenergy_use_rule)

电解器始终受到电解剂输入功率的限制

def H2energy_production_rule(model, l, h):
    return model.E_h2[l, h] <= model.Ph2[l]
model.H2energy_production = Constraint(model.L, model.H, rule=H2energy_production_rule)

氢收率等于每个位置的效率

def H2hydrogen_production_rule(model, l, h):
    return model.Y_h2[l, h] == model.E_h2[l, h] * effi
model.H2hydrogen_production = Constraint(model.L, model.H, rule=H2hydrogen_production_rule)

在每个位置，两者的总安装能力是发电厂和电解器不得超过可用区域的

def power_density_rule(model,l):
    return model.Psolar[l]/Pden_solar + model.Pwind[l]/Pden_wind + model.Ph2[l]/Pden_h2 <= alpha
model.power_density = Constraint(model.L, rule=power_density_rule)

目标函数

最小化设置中所有点

def objective_function(model):
    return sum((K_solar*model.Psolar[l] + K_wind*model.Pwind[l] + K_h2*model.Ph2[l]) / ((sum(model.Y_h2[l,h] for h in model.H))*Sy) for l in model.L)
model.obj = Objective(rule=objective_function, sense=minimize)

求解模型

求解的总和，但应与任何NLP求解器一起使用（例如IPOPT）

instance = model.create_instance(data)
#optimizer=SolverFactory('ipopt')
optimizer = SolverFactory('octeract-engine', executable=r"C:\Program Files\Octeract\bin\octeract-engine.exe")
results = optimizer.solve(instance, tee = True)

获得结果的状态

status = results.solver.status
print("Status:", status)
termination = results.solver.termination_condition
print("Termination:", termination)

model.pprint()

一个位置的错误消息

，求解器运行约1分钟（对于五个位置，几乎是10分钟），然后说明：

退出：超过最大CPU时间。
警告：将带有警告状态的求解对象加载到模型=未知;
状态：警告
终止：Maxiterations

对于这两个测试，model.pprint（）也输出模型的内容，这些内容始终均以 size = 0 and 和未构建<< /em>。

导出数据

将决策变量导出到CSV；每次结果为零。

Psolar=instance.Psolar
Psolar_out = pd.DataFrame(columns=[['id1','id2']])
for index in Psolar:
    Psolar_out.at[index,'id1']=value(Psolar[index])
Psolar_out.to_csv('_psolar.csv')


Pwind=instance.Psolar
Pwind_out = pd.DataFrame(columns=[['id1','id2']])
for index in Pwind:
    Pwind_out.at[index,'id1']=value(Pwind[index])
Pwind_out.to_csv('_pwind.csv')


Ph2=instance.Psolar
Ph2_out = pd.DataFrame(columns=[['id1','id2']])
for index in Ph2:
    Ph2_out.at[index,'id1']=value(Ph2[index])
Ph2_out.to_csv('_ph2.csv')

Answer 1

可能会有所帮助的几个想法/观察结果：

1. 您的“减少功率”的变量可能会被删除，只会使相关的约束不平等。那将消除一些看似不必要的（？）变量。例如，如果是：

Tot_power = Hydrogen_power + curtailed_power

to：

Tot_power >= Hydrogen Power

2. 对于非线性模型，您的模型非常大。我只是戳了数据。在非线性模型中解决（什么是什么？）10,000个时间段可能太多了。您可以：使您的模型线性（好主意）或每天或6小时的块汇总，或者在不同时期内运行模型并比较答案。

3. 使您的模型线性启动并通过删除该部门来使其正常工作。看看那是你的地方。线性模型应几乎立即解决。

4. 模型中似乎没有任何激励/压力来制造任何力量。您正在限制上层生产，但是您没有我能看到的需求信号，并且您正在最小化成本，因此您可能会得到默认解决方案：毫无疑问/什么都不做==＆gt;成本为零。 ：）

5. 仅在求解器中获得良好的状态，仅检查数据值。任何其他错误/不可行/超时/等。意味着加载的数据是垃圾。

6. 您可以通过将目标从比率（使用变量分配）更改为对除数的惩罚来实验并使模型线性：

   sum（p_wind + p_solar + p_h2） - w * sum（y_h2）

其中是一个合理的权衡价值。 +某些需求信号（或利润）应获得模型呼吸。记住：

linear_model_that_works_but_has_some_warts >> broken_model

Answer 2

从快速查看您的数据来看，大小的尺寸就会差不多。 Octeract Engine具有最先进的前置主义器，因此最好将其放在发动机上以消除变量，因为它将以数值稳定的方式进行。如果引擎确定可以用来收紧问题的放松，则该引擎还可以决定保留一些冗余的约束和变量，因此最好将重点放在模型中的逻辑上，而不必担心冗余信息。

当求解器在预设阶段及其分支和绑定之前，我已经看到了特定的PYOMO消息。当求解器不产生.sol file时，这是一个pyomo的限制，但这是实际发生的事情 - 没有求解器日志，很难说，但是我假设您已经指定了某些终止条件的条件求职者超时。

如果没有，唯一想到的是引擎将您的问题视为凸面，将其传递给iPopt，然后iPopt无法解决。在这种情况下，这是预期的行为，因为发动机本身可以做到。

您可以尝试尝试的一个选项是指定use_approx_equaltial_constraints = true在OCTERACT-INGINE选项中，这将重新重新重新制定您的平等约束。虽然这在一般情况下是一个坏主意，但我添加了该选项，因为它通常有助于ipopt在可行性中挣扎的大型NLP。

您还可以尝试更改起点选择策略。

此外，您可以通过在执行目录中添加ipopt.opt.opt.opt.opt选项文件来控制OCTERACT-INGINE所调用的IPOPT。这将使您能够准确控制IPOPT正在做的事情，包括更改其约束违规选项。