问如何在Pyomo中对两阶段随机规划中的情景进行适当的索引？

EN

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`