基于Xpress库的船舶Retrofit优化约束问题技术求助

问题背景 #

将Excel文件数据读入DataFrame，建立优化模型最小化16艘船舶在2023-2050年（共27年）的总排放。船舶改造后可永久减排5%，需确定每艘船的改造时间（或是否改造）以最小化总排放。当前约束逻辑无法实现“改造后永久享受减排收益”的要求，需修正。

原问题约束代码段 #

# Constrain the binary variables is_retrofitted and retrofitted so that they are consistent
for i in vessels:
    m.addConstraint(retrofitted[i] <= is_retrofitted[np.ravel_multi_index((i, year_retrofitted[i]), (len(vessels), len(time)))])
    m.addConstraint(xp.Sum(is_retrofitted[i, t] for t in time) == retrofitted[i])

原完整代码 #

# Importing packages
import xpress as xp
import pandas as pd
import numpy as np

# Create path to excel file with fleet information
path = r'path_to_excel_file'
df = pd.read_excel(path, sheet_name = 'sheet_name') # Create dataframe from excel sheet
df.index = df['Vessel'] # Set index equal to the vessels

# Adjusting the dataframe
df = df.drop('Total') # Delete unnecessary row
df = df.drop('Unit') # Delete unnecessary row
df.drop('Total', axis=1, inplace=True) # Delete unnecessary column

# Define sets
V = df['Vessel'] # Set of vessels
df = df.drop('Vessel', axis = 1) # Delete unnecessary column

# Define the number of vessels and the time horizon
num_vessels = 16
time_horizon = 2050 - 2023 + 1

# Create sets for vessels and time
vessels = range(num_vessels)
time = range(time_horizon)

# Create the emissions parameter
emissions = [xp.var("emissions_{}_{}".format(i, t), lb=0) for i in vessels for t in time]

# Create indices to map each value to a specific vessel and year
index = [(i, t) for i in vessels for t in time]

# Iterate over the rows and columns of the dataframe
for (i, t), e in zip(index, emissions):
    e.lb = df.iloc[i, t + 2023 - 2023]

# Define a binary variable for each vessel, indicating whether or not it is retrofitted
retrofitted = [xp.var("retrofitted_{}".format(i), vartype=xp.binary) for i in vessels]

# Define a binary variable for each vessel and year, indicating whether or not the vessel is retrofitted in that year
is_retrofitted = [xp.var("is_retrofitted_{}_{}".format(i, t), vartype=xp.binary) for i in vessels for t in time]

# Define an integer variable for each vessel, indicating the year in which it is retrofitted
year_retrofitted = [xp.var("year_retrofitted_{}".format(i), vartype=xp.integer) for i in vessels]

# Create the objective function to minimize the total CO2 emissions of the vessel fleet
m = xp.problem("Minimize CO2 emissions")

# Add variable to problem
m.addVariable(emissions)
m.addVariable(retrofitted)
m.addVariable(is_retrofitted)
m.addVariable(year_retrofitted)

# Set objective funtion
m.setObjective(xp.Sum(emissions[np.ravel_multi_index((i, t), (len(vessels), len(time)))] * (1 - 0.05 * is_retrofitted[i])
                      for i, vessel in enumerate(vessels) for t, year in enumerate(time)), sense=xp.minimize)

# Constrain the binary variables is_retrofitted and retrofitted so that they are consistent
for i in vessels:
    m.addConstraint(retrofitted[i] <= is_retrofitted[np.ravel_multi_index((i, year_retrofitted[i]), (len(vessels), len(time)))])
    m.addConstraint(xp.Sum(is_retrofitted[i, t] for t in time) == retrofitted[i])

# Constrain the year_retrofitted variables to be between 0 (indicating no retrofit) and the time horizon
for i in vessels:
    m.addConstraint(year_retrofitted[i] >= 0)
    m.addConstraint(year_retrofitted[i] <= time_horizon)

# Solve the optimization problem
m.solve()

# Get the optimal values of the binary variables
is_retrofitted_opt = [is_retrofitted[i].getValues() for i in vessels]
retrofitted_opt = [retrofitted[i].getValues() for i in vessels]
year_retrofitted_opt = [year_retrofitted[i].getValues() for i in vessels]

# Get the optimal value of the objective function
obj_value = m.getObjVal()

# Print the results
print("The minimum total CO2 emissions of the vessel fleet is:", obj_value)
print("The optimal retrofit decisions are:", retrofitted_opt)
print("The optimal years for retrofitting are:", year_retrofitted_opt)

问题分析 #

原代码存在3个核心问题：

1. 变量使用错误：目标函数中引用is_retrofitted[i]（缺少时间维度），无法区分船舶在不同年份是否已完成改造
2. 索引逻辑错误：尝试用变量year_retrofitted[i]作为np.ravel_multi_index的参数，这在优化建模中不合法（索引必须是常量）
3. 约束逻辑缺失：未设置“一旦改造，后续所有年份持续享受减排”的约束，无法实现永久减排的业务规则

修正方案与代码 #

重新定义变量并补充关键约束，实现改造后永久生效的逻辑：

修正后的完整代码 #

# Importing packages
import xpress as xp
import pandas as pd
import numpy as np

# Create path to excel file with fleet information
path = r'path_to_excel_file'
df = pd.read_excel(path, sheet_name='sheet_name')
df.index = df['Vessel']

# Adjusting the dataframe
df = df.drop('Total')
df = df.drop('Unit')
df.drop('Total', axis=1, inplace=True)

# Define sets
V = df['Vessel']
df = df.drop('Vessel', axis=1)

# Define the number of vessels and the time horizon
num_vessels = 16
time_horizon = 2050 - 2023 + 1
vessels = range(num_vessels)
time = range(time_horizon)

# Create emissions parameter (base emissions without retrofit)
base_emissions = df.values  # Directly get base emission values from dataframe

# --- Variable Definitions ---
# Binary: 1 if vessel i is ever retrofitted
retrofitted = [xp.var(f"retrofitted_{i}", vartype=xp.binary) for i in vessels]
# Binary: 1 if vessel i is retrofitted *in year t*
is_retrofitted_in_year = [[xp.var(f"is_retrofitted_in_year_{i}_{t}", vartype=xp.binary) 
                           for t in time] for i in vessels]
# Integer: Year when vessel i is retrofitted (0 to time_horizon; set to time_horizon if not retrofitted)
retrofit_year = [xp.var(f"retrofit_year_{i}", vartype=xp.integer, lb=0, ub=time_horizon) for i in vessels]
# Binary: 1 if vessel i has been retrofitted by year t (enjoys 5% emission reduction)
is_retrofit_active = [[xp.var(f"is_retrofit_active_{i}_{t}", vartype=xp.binary) 
                       for t in time] for i in vessels]

# Create optimization problem
m = xp.problem("Minimize CO2 emissions")
m.addVariable(retrofitted, is_retrofitted_in_year, retrofit_year, is_retrofit_active)

# --- Objective Function ---
# Minimize total emissions: base_emissions[i,t] * (1 - 0.05 * is_retrofit_active[i,t])
total_emissions = xp.Sum(base_emissions[i][t] * (1 - 0.05 * is_retrofit_active[i][t])
                         for i in vessels for t in time)
m.setObjective(total_emissions, sense=xp.minimize)

# --- Constraints ---
# 1. Each vessel can be retrofitted at most once (0 or 1 times)
for i in vessels:
    m.addConstraint(xp.Sum(is_retrofitted_in_year[i][t] for t in time) == retrofitted[i])

# 2. Link retrofit_year to is_retrofitted_in_year: if retrofitted, year is the t where is_retrofitted_in_year[i][t]=1
#    If not retrofitted, set retrofit_year to time_horizon (outside our time range)
for i in vessels:
    m.addConstraint(retrofit_year[i] == xp.Sum(t * is_retrofitted_in_year[i][t] for t in time) + 
                    (1 - retrofitted[i]) * time_horizon)

# 3. Ensure is_retrofit_active[i][t] = 1 if retrofit_year[i] <= t (retrofitted by year t)
for i in vessels:
    for t in time:
        # If is_retrofit_active[i][t] = 1, then retrofit_year[i] <= t
        m.addConstraint(retrofit_year[i] <= t + time_horizon * (1 - is_retrofit_active[i][t]))
        # If not retrofitted, is_retrofit_active must be 0
        m.addConstraint(is_retrofit_active[i][t] <= retrofitted[i])

# 4. Recursive constraint: once active, stay active (optional but reinforces logic)
for i in vessels:
    # Initial year: active if retrofitted in year 0
    m.addConstraint(is_retrofit_active[i][0] == is_retrofitted_in_year[i][0])
    # For subsequent years: active if previous was active OR retrofitted this year
    for t in range(1, time_horizon):
        m.addConstraint(is_retrofit_active[i][t] >= is_retrofit_active[i][t-1])
        m.addConstraint(is_retrofit_active[i][t] >= is_retrofitted_in_year[i][t])
        m.addConstraint(is_retrofit_active[i][t] <= is_retrofit_active[i][t-1] + is_retrofitted_in_year[i][t])

# --- Solve and Output ---
m.solve()

# Extract results
retrofitted_opt = [var.getValues() for var in retrofitted]
retrofit_year_opt = [var.getValues() for var in retrofit_year]
is_retrofit_active_opt = [[var.getValues() for var in row] for row in is_retrofit_active]

# Calculate total emissions
obj_value = m.getObjVal()

# Print results
print(f"最小船队总CO2排放量: {obj_value}")
print(f"最优改造决策 (1=改造, 0=不改造): {retrofitted_opt}")
print(f"最优改造年份 (0=2023, ..., 27=2050; 27表示不改造): {retrofit_year_opt}")

关键修正说明 #

1. 变量拆分：将原is_retrofitted拆分为is_retrofitted_in_year（当年改造）和is_retrofit_active（当年已改造），明确区分两种状态
2. 合法索引关联：通过线性约束关联retrofit_year和is_retrofitted_in_year，避免用变量作为索引的错误
3. 永久生效约束：通过retrofit_year[i] <= t + time_horizon*(1 - is_retrofit_active[i][t])确保改造后所有年份is_retrofit_active为1，实现永久减排
4. 目标函数修正：使用is_retrofit_active[i][t]计算当年减排量，逻辑准确

内容的提问来源于stack exchange，提问作者Felix Dietrichson