如何在 Pyomo 中的两阶段随机编程中正确索引场景？

问题描述 投票：0回答：1

主要问题： 如何定义风场景参数，使其能够针对场景集“S”中不同的概率提供不同的风场景。

运行场景集 S = [1,2,3,...10] 且相应概率 Prob = [0.1, 0.1, ..., 0.1] 的随机代码时。 t = [1, 2, ..., 48] 以及相应的风速 m.wind.

我尝试了以下选项来为下面的目标函数创建 (s*t) 矩阵。

Objective function

选项1： m.wind = dict{:48}，这是确定性模型中使用的 48 个时间步长的风速字典

选项2： 米。 Wind = {ndarray:(10,48)} = array([[...],..,[...]]), 在此我构造了一个包含 10 个场景的数组，其中包括 10 个时间步长的风速场景。

选项3： m.windd = {list:10} = array[(...),...,(...)], 在此选项中，在读到 Pyomo 有时无法识别正方形后，我将不同的场景放在括号之间括号。

选项4： 构建风场景集的最后一种方法是创建 (x*t) 字典。

所有选项都导致错误： 错误：索引“0”对于组件“wind”确实无效

您知道如何解决此错误以及如何正确索引风吗？

`

def build_model(price_data, horizon_length, scenario_length, load_calc, park_calc):

m = pyo.ConcreteModel()

### BEGIN SOLUTION

# test vector
vector = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1])
vector = vector.reshape(10,1)

## Sets

# Save the number of timesteps
m.N = horizon_length
m.S = len(scenario_length)

# Define the horizon set starting at hour 1 until horizon length +1
m.HORIZON = pyo.Set(initialize=range(1, m.N + 1))

# Define scenario set
m.SCENARIO = pyo.Set(initialize=range(1, m.S + 1))

## Parameters

# Round trip efficiency
m.teta = pyo.Param(initialize=0.95)

# Energy [MWh] in battery at t=0
m.E0 = pyo.Param(initialize=2.0, mutable=True)

# Guarantee of origin for local wind [€/MWh]
m.goNL = pyo.Param(initialize=5)

# Guarantee of origin for grid power [€/MWh]
m.goBE = pyo.Param(initialize=150)

# Maximum discharge power
m.d_max = pyo.Param(initialize=5)

# Maximum charge power
m.c_max = pyo.Param(initialize=5)

# Maximum export power
m.im_max = pyo.Param(initialize=10)

# Maximum import power
m.ex_max = pyo.Param(initialize=100)

## CREATE DICTS FOR DATA: Price, Load & Calc
# Create empty dictionary
price_data_dict = {}

# Loop over Price data elements of numpy array
for i in range(0, len(price_data)):
    # Add element to data_dict
    price_data_dict[i + 1] = price_data[i]

# Create empty dictionary
load_data_dict = {}

# Loop over Load data elements of numpy array
for i in range(0, len(load_calc)):
    # Add element to data_dict
    load_data_dict[i + 1] = load_calc[i]

# Create empty dictionary
park_data_dict = {}

# Loop over Wind park data elements of numpy array
for i in range(0, len(park_calc)):
    # Add element to data_dict
    park_data_dict[i + 1] = park_calc[i]

# Create empty dictionary
prob_dict = {}

# Loop over probability data elements of numpy array
for i in range(0, len(vector)):
    # Add element to prob_dict
    prob_dict[i + 1] = vector[i]

# Repeat the wind data to a matrix for 10 similar scenario's
wind_matrix = np.tile(park_calc, (10, 1))
# wind_matrix = np.tile(park_calc, (10, 1)) * vector

park_data_dict_2 = {1: park_data_dict, 2: park_data_dict, 3: park_data_dict, 4: park_data_dict, 5: park_data_dict,
                    6: park_data_dict, 7: park_data_dict, 8: park_data_dict, 9: park_data_dict, 10: park_data_dict}

# Price data
m.price = pyo.Param(m.HORIZON, initialize=price_data_dict, domain=pyo.Reals, mutable=True)

# Load data
m.Load = pyo.Param(m.HORIZON, initialize=load_data_dict, domain=pyo.Reals, mutable=True)

# Wind park data
m.wind = pyo.Param(m.SCENARIO, m.HORIZON, initialize=park_data_dict_2, mutable=True) #park_data_dict

# Scenario probability
m.prob = pyo.Param(m.SCENARIO, initialize=vector)  # Was Scen_prob

# # New description of wind in 10 different scenarios
# m.wind = pyo.Param(m.SCENARIO, m.HORIZON, initialize=wind_matrix_2) #  initialize=wind_matrix_2

## Variables

## Battery related variables
# Charging rate [MW]
m.c = pyo.Var(m.HORIZON, initialize=0.0, bounds=(0, 10), domain=pyo.NonNegativeReals)

# Discharging rate [MW]
m.d = pyo.Var(m.HORIZON, initialize=0.0, bounds=(0, 10), domain=pyo.NonNegativeReals)

# Battery power
m.Bat = pyo.Var(m.HORIZON, initialize=0.0, domain=pyo.NonNegativeReals)

# Binary variables charging and grid
m.u = pyo.Var(m.HORIZON, initialize=0.0, domain=pyo.Binary)
m.v = pyo.Var(m.HORIZON, initialize=0.0, domain=pyo.Binary)

# Energy (state-of-charge) [MWh]
m.E = pyo.Var(m.HORIZON, initialize=2.0, bounds=(0, 5), domain=pyo.NonNegativeReals)
m.G_im = pyo.Var(m.HORIZON, initialize=0, bounds=(0, 10), domain=pyo.NonNegativeReals)
m.G_ex = pyo.Var(m.HORIZON, initialize=0, bounds=(0, 100), domain=pyo.NonNegativeReals)
m.grid = pyo.Var(m.HORIZON, initialize=m.Load, bounds=(0, 10), domain=pyo.NonNegativeReals)

# Objective function

# def objfun(model):
#     return sum((m.price[t] + m.goNL) * m.wind[t] + (m.price[t] + m.goBE) * m.G_im[t] for t in m.HORIZON)

def objfun(model):
    return sum((m.price[t] + m.goBE) * m.G_im[t] + (m.price[t] + m.goNL) * sum(m.prob[s] * m.wind[s, t] for s in m.SCENARIO) for t in m.HORIZON)

m.OBJ = pyo.Objective(rule=objfun, sense=pyo.minimize)

def PowerBalance(m, t):
    return m.Load[t] + m.c[t] == m.grid[t] + m.d[t]

# Define Energy Balance constraints. [MWh] = [MW]*[1 hr]
# Note: assume 1-hour timestep in price data and control actions.
def EnergyBalance(m, t):
    # First timestep
    if t == 1:
        return m.E[t] == m.E0 + m.c[t] * m.teta - m.d[t] / m.teta

        # Subsequent timesteps
    else:
        return m.E[t] == m.E[t - 1] + m.c[t] * m.teta - m.d[t] / m.teta

# def ColdIroning(m, t):
#     return m.c[t] + m.d[t] + m.Load[t] <= m.CI

def GridBalance(m, t, s):
    return m.grid[t] == m.wind[t, s] + m.G_im[t] - m.G_ex[t]

def ImMax(m, t):
    return m.G_ex[t] - m.v[t] * m.ex_max <= 0

def ExMax(m, t):
    return m.G_im[t] + m.v[t] * m.im_max <= m.im_max

# def BatteryBalance(m, t):
#     return m.Bat[t] - m.d[t] + m.c[t] == 0
#
def ChargeMax(m, t):
    return m.d[t] - m.u[t] * m.d_max <= 0

def DischargeMax(m, t):
    return m.c[t] + m.u[t] * m.c_max <= m.c_max

m.EnergyBalance_Con = pyo.Constraint(m.HORIZON, rule=EnergyBalance)
m.PowerBalance_Con = pyo.Constraint(m.HORIZON, rule=PowerBalance)
# m.ColdIroning_Con = pyo.Constraint(m.HORIZON, rule=ColdIroning)
m.GridBalance_Con = pyo.Constraint(m.HORIZON, m.SCENARIO, rule=GridBalance)
# m.BatteryBalance_Con = pyo.Constraint(m.HORIZON, rule=BatteryBalance)
m.ChargeMax_Con = pyo.Constraint(m.HORIZON, rule=ChargeMax)
m.DischargeMax_Con = pyo.Constraint(m.HORIZON, rule=DischargeMax)
m.ImMax_Con = pyo.Constraint(m.HORIZON, rule=ImMax)
m.ExMax_Con = pyo.Constraint(m.HORIZON, rule=ExMax)
## END SOLUTION

return m`
python pyomo dispatch scenarios stochastic
1个回答
0
投票

模型上的两件事：

您不需要（也可能不应该）初始化变量，只需让求解器完成其工作即可。

你的 

GridBalance
约束中有一个严重的拼写错误 Wind[t, s] （而不是 Wind[s, t]），如果 |T| 找到的话，这将是一个魔鬼。 == |S|。

你没有说风数据以什么格式传给你。也许你只是手工插入一些表格数据。那么让我们从 

pyomo
想要什么开始吧。它需要一个元组索引字典来初始化参数。这意味着字典的键值是 (s, t) 值的元组。这也称为“平面”数据结构，其中所有键都被枚举，并且感兴趣的数据值位于 1 列中（即矩阵格式或其他格式）。

所以你想用这样的方式初始化你的参数：

import pyomo.environ as pyo

# what we want: a "flat" dictionary

#         s, t : w
wind = { (1, 1): 12,
         (1, 2): 11,
         (1, 3): 10,
         (2, 1): 9,
         (2, 2): 13,
         (2, 3): 14}

m = pyo.ConcreteModel()

# SETS
m.S = pyo.Set(initialize=[1, 2])
m.T = pyo.Set(initialize=[1, 2, 3])

# PARAMS

m.wind = pyo.Param(m.S, m.T, initialize=wind)

m.pprint()

输出：

3 Set Declarations
    S : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {1, 2}
    T : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    3 : {1, 2, 3}
    wind_index : Size=1, Index=None, Ordered=True
        Key  : Dimen : Domain : Size : Members
        None :     2 :    S*T :    6 : {(1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3)}

1 Param Declarations
    wind : Size=6, Index=wind_index, Domain=Any, Default=None, Mutable=False
        Key    : Value
        (1, 1) :    12
        (1, 2) :    11
        (1, 3) :    10
        (2, 1) :     9
        (2, 2) :    13
        (2, 3) :    14

4 Declarations: S T wind_index wind

显然，有多种技术可以根据数据创建元组索引字典，具体取决于源数据的结构。但在此之前，我看到您使用循环将列表转换为字典。当然可行，或者您可以使用 enumerate 的快捷方式，它生成索引：值对并将其传递给字典构造函数。请注意，如果您喜欢数据 1 索引副 0，则 enumerate 采用可选的起始参数。

prices = [3.5, 4.2, 9.8]
price_dict = dict(enumerate(prices, start=1))

print(price_dict)

# {1: 3.5, 2: 4.2, 3: 9.8}

因此，如果您纯粹手动输入值，您可以输入如上所示的平面字典，或者如果您有风数据的列表列表（又称矩阵），您可以通过多种方式对其进行转换，具体取决于您对字典理解等的舒适程度。以下所有 3 个都会生成可在您的模型中使用的相同平面字典：

raw_wind_data = [[4, 5, 9],
                 [3, 0, 12]]

wind_1 = {}
for i in range(len(raw_wind_data)):
    for j in range(len(raw_wind_data[0])):
        wind_1[(i+1, j+1)] = raw_wind_data[i][j]

wind_2 = { (r+1, c+1) : raw_wind_data[r][c] 
                    for r in range(len(raw_wind_data)) 
                    for c in range(len(raw_wind_data[0]))}

wind_3 = {(r_idx, c_idx): w 
            for r_idx, row in enumerate(raw_wind_data, 1)
            for c_idx, w   in enumerate(row, 1)}

print(wind_1)
print(wind_2)
print(wind_3)

# {(1, 1): 4, (1, 2): 5, (1, 3): 9, (2, 1): 3, (2, 2): 0, (2, 3): 12}
# {(1, 1): 4, (1, 2): 5, (1, 3): 9, (2, 1): 3, (2, 2): 0, (2, 3): 12}
# {(1, 1): 4, (1, 2): 5, (1, 3): 9, (2, 1): 3, (2, 2): 0, (2, 3): 12}
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.