将优化器的目标设置为输入的标准偏差（使用PYMO的非线性优化）

Question

我正在尝试将PYMO用于单个目标非线性优化问题。

目标函数是最大程度地减少输入变量的方差（或标准偏差）之后的某些约束（我可以在Excel中进行）。

以下是我想做的错误的代码

model = pyo.ConcreteModel()

# declare decision variables
model.x1 = pyo.Var(domain=pyo.NonNegativeReals)
model.x2 = pyo.Var(domain=pyo.NonNegativeReals)
model.x3 = pyo.Var(domain=pyo.NonNegativeReals)
model.x4 = pyo.Var(domain=pyo.NonNegativeReals)

# declare objective
from statistics import stdev
model.variance = pyo.Objective(
    expr = stdev([model.x1, model.x2, model.x3, model.x4]),
    sense = pyo.minimize)

# declare constraints
model.max_charging = pyo.Constraint(expr = model.x1 + model.x2 + model.x3 + model.x4 >= 500)
model.max_x1 = pyo.Constraint(expr = model.x1 <= 300)
model.max_x2 = pyo.Constraint(expr = model.x2 <= 200)
model.max_x3 = pyo.Constraint(expr = model.x3 <= 100)
model.max_x4 = pyo.Constraint(expr = model.x4 <= 200)

# solve
pyo.SolverFactory('glpk').solve(model).write()

#print
print("energy_price = ", model.variance())
print(f'Variables = [{model.x1()},{model.x2()},{model.x3()},{model.x4()}]')

。

示例 stdev 功能来自统计信息。

我的假设是模型变量x1-x4尚未分配一个值，这是主要问题。但是，我不确定如何处理？

Answer 1

第一：stdev是非线性的。那么，为什么还要尝试用线性求解器解决这个问题呢？

PYOMO不知道统计信息包。您必须使用基本操作，使用外部功能方法或使用近似值（例如最小化范围）来对标准偏差进行编程。

Answer 2

因此，我设法解决了这个问题，并且包括以下解决方案。但是首先，正如@erwin Kalvelagen提到的那样，我想突出一些要点

1. ，“ stdev”是非线性的，因此像'glpk'这样的线性求解器总是会导致错误。就我的问题而言，“ iPopt”工作正常，但要小心，因为在某些情况下，它的表现会很差。
2. 另外，正如@erwin Kalvelagen所提到的那样，Pyomo对统计软件包不了解。因此，当您尝试使用该软件包（例如'stdev'，“方差”等）中的函数时，它将尝试在求解器分配任何值之前尝试评估模型变量并导致错误。
3. PYOMO目标函数需要表达。下面的代码示例显示了如何在不使用外部软件包的情况下动态生成方差的表达式。代码对您拥有的模型变量的数量也不可知。
4. 使用方差或标准偏差将为我的项目提供相同的目的。我选择使用方差来避免计算其平方根（因为标准偏差是方差的平方根）。

可变性目标函数

import pyomo.environ as pyo
def variabilityRule(model):
    #Calculate mean of all model variables
    for index in model.x:
        mean += model.x[index]
    
    mean = mean/len(model.x)

    #Calculate the difference between each model variable and the mean
    obj_exp = ((model.x[1])-mean) * ((model.x[1])-mean)
    for i in range(2, len(model.x)+1): #note that pyomo variables start from 1 not zero
        obj_exp += ((model.x[i])-mean) * ((model.x[i])-mean)

    #Divide by number of items
    obj_exp = (obj_exp/(len(model.x)))

    return obj_exp

model = pyo.ConcreteModel()
model.objective = pyo.Objective(
    rule = variabilityRule,
    sense = pyo.maximize)

编辑：标准偏差计算

You can use the same approach to calculate the standard deviation as I found out. Just multiply the final expression (`obj_exp`) by power 0.5

obj_exp = obj_exp ** 0.5

...

ps如果您有兴趣，这就是我基于输入数组动态生成模型变量的方式

model.x = pyo.VarList(domain=pyo.NonNegativeReals)
for i in range(0, len(input_array)):
    model.x.add()