基于Gekko的光伏EMS系统MINLP优化：产量约束下利润最大化

我正在搭建一套光伏系统的能源管理系统（EMS）原型，最初用MILP求解器，但因目标函数构建导致问题转为MINLP，改用Gekko后表现不错。

当前系统需根据电解制氢情况，决策是否购电生产绿氨——电解制氢和哈伯法合成氨的能源来自光伏或电网（电网电需购买），同时要决策氨的销售市场以最大化利润。

现在要加约束：月度氨产量至少30吨（30000kg），仅靠光伏无法满足，必须购电才能达标，但同时还要最大化利润。

我试过用双目标函数分别最大化产量和利润，但Gekko优先考虑利润，根本不购电。我的问题是：能不能让Gekko先达成指定产量，再最优实现利润最大化？

我考虑过用条件判断切换目标函数：先最大化产量，达标后再切换到利润最大化，但这似乎不是最优方案。另外，我用while循环逐小时做月度优化，不确定这是不是Gekko里最规范的实现方式。

以下是我目前的代码（已加入条件判断逻辑）：

#Attempt using Gekko
from gekko import GEKKO
i=0 #counter for the while loop
contract_prod=30000 #kg

nh3_ue=[]#measures the percentage of ammonia sold to the European Union
nh3_usa=[]#measures the percentage of ammonia sold to the US
nh3_asia=[]#measures the percentage of ammonia sold to ASIA
nh3_br=[]#measures the percentage of ammonia sold to Brazil
max_profit=[]#measures the profit
p_buy = []#measures how much energy was bought
accumulate_prod=[]#measures accumulated production
prod_index=[]#measures the production in a specific hour

while i<=24:
#Initialize Model
  m = GEKKO()

#Set Global Options
  m.options.SOLVER=1 

# optional solver settings with APOPT
  m.solver_options = ['minlp_maximum_iterations 500', \
                    # minlp iterations with integer solution
                    'minlp_max_iter_with_int_sol 10', \
                    # treat minlp as nlp
                    'minlp_as_nlp 0', \
                    # nlp sub-problem max iterations
                    'nlp_maximum_iterations 50', \
                    # 1 = depth first, 2 = breadth first
                    'minlp_branch_method 1', \
                    # maximum deviation from whole number
                    'minlp_integer_tol 0.05', \
                    # covergence tolerance
                    'minlp_gap_tol 0.01']

#Parameters
  c_H2 = m.Param(value=4.5,name="LCOH")
  c_HB = m.Param(value=0.757,name="xCost_Haber-Bosch")
  p_ue = m.Param(value=1.602,name="Price_europe")
  p_usa = m.Param(value=1.323,name="Price_usa")
  p_asia = m.Param(value=0.84,name="Price_asia")
  p_br=m.Param(value=1.55,name="Price_brazil")
  per_min_br=m.Param(value=20,name="minimmum_brasil")
  perc=m.Param(value=100,name="total_percentage_ammonia")
  limit=m.Param(value=963070.78)#Electrolyzer limit
  p_ufvi=m.Param(p_ufv[0][i])#Energy from PV in ith hour
  c_Ci=m.Param(c_C[0][i])#Grid energ cost in ith hour

#Initialize Variables
  per_eur = m.Var(value=2,lb=0,ub=100,name="ammonia_to_europe_percentage")
  per_usa = m.Var(value=2,lb=0,ub=100,name="ammonia_to_usa_percentage")
  per_asia = m.Var(value=2,lb=0,ub=100,name="ammonia_to_asia_percentage")
  per_br= m.Var(value=2,lb=0,ub=100,name="ammonia_to_br_percentage")
  p_C = m.Var(value=100,lb=0,ub=963070.78,name="energy_bought")

#Equations
  m.Equation(per_eur+per_usa+per_asia+per_br+per_min_br<=perc)
  m.Equation(p_ufvi+p_C<=limit)

#Objective
  obj1 = m.Intermediate(per_eur*(p_ue-c_HB-c_H2/5.67)*(p_ufvi+p_C)/1108500
                        +per_usa*(p_usa-c_HB-c_H2/5.67)*(p_ufvi+p_C)/1108500
                        +per_asia*(p_asia-c_HB-c_H2/5.67)*(p_ufvi+p_C)/1108500
                        +(per_br+per_min_br)*(p_br-c_HB-c_H2/5.67)*(p_ufvi+p_C)/1108500
                        -c_Ci*p_C/1000000) #Profit
  obj2=m.Intermediate((p_ufvi+p_C)/11085) #Production

#Objective
  if i==0:
    m.Maximize(obj2)
  elif (accumulate_prod[i-1]<=contract_prod):
    m.Maximize(obj2)
  else:#maximizar lucro
    m.Maximize(obj1)
 
#Open the folder created with the results
  m.open_folder()

#Solve simulation
  try:
    m.solve(disp=True)    # solve
  except:
    print('Not successful')
    from gekko.apm import get_file 
    print(m._server)
    print(m._model_name)
    f = get_file(m._server,m._model_name,'infeasibilities.txt')
    f = f.decode().replace('\r','')
    with open('infeasibilities.txt', 'w') as fl:
        fl.write(str(f))

#Results
  nh3_ue.append(per_eur.value[0])
  nh3_usa.append(per_usa.value[0])
  nh3_asia.append(per_asia.value[0])
  nh3_br.append(per_br.value[0])
  p_buy.append(p_C.value[0])
  max_profit.append(obj1.value[0])
  prod_index.append(obj2.value[0])
  if i==0:
    acc=obj2.value[0]
  else:
    j=i-1
    acc=accumulate_prod[j]+obj2.value[0]
  accumulate_prod.append(acc)
  i=i+1

补充参数（24小时数据） #

p_ufv（光伏发电量，单位W）

[0        0.0
 1        0.0
 2        0.0
 3        0.0
 4        0.0
 5        0.0
 6        0.0
 7        78843.24
 8        330970.17
 9        558981.33
 10       735160.91
 11       800000.00
 12       800000.00
 13       800000.00
 14       800000.00
 15       755987.50
 16       587748.41
 17       366223.30
 18       113380.86
 19       0.0
 20       0.0
 21       0.0
 22       0.0
 23       0.0
Name: Potência W, Length: 24, dtype: float64]

c_C（电网购电成本）

[0        379.39
 1        379.39
 2        379.39
 3        379.39
 4        379.39
 5        379.39
 6        379.39
 7        379.39
 8        379.39
 9        379.39
 10       379.39
 11       379.39
 12       379.39
 13       379.39
 14       379.39
 15       379.39
 16       379.39
 17       527.01
 18       527.01
 19       527.01
 20       527.01
 21       379.39
 22       379.39
 23       379.39
Name: Preco, Length: 24, dtype: float64]

1. 优先满足产量约束的正确做法 #

不要用逐小时切换目标函数的方式，这会导致局部最优——每一步只看当前小时的决策，没有全局统筹。正确的做法是把产量要求作为硬约束加入模型，直接最大化利润，求解器会自动在满足产量要求的前提下，找到利润最高的方案。

2. 重构为全局优化模型（替代逐小时循环） #

逐小时循环的核心问题是决策独立，无法利用跨时段的成本差异（比如在电价低时多购电生产，既满足配额又降低成本）。重构后的全局模型示例如下：

from gekko import GEKKO
import pandas as pd

# 加载参数数据
p_ufv = pd.Series([0.0,0.0,0.0,0.0,0.0,0.0,0.0,78843.24,330970.17,558981.33,
                   735160.91,800000.00,800000.00,800000.00,800000.00,755987.50,
                   587748.41,366223.30,113380.86,0.0,0.0,0.0,0.0,0.0])
c_C = pd.Series([379.39,379.39,379.39,379.39,379.39,379.39,379.39,379.39,
                 379.39,379.39,379.39,379.39,379.39,379.39,379.39,379.39,
                 379.39,527.01,527.01,527.01,527.01,379.39,379.39,379.39])
contract_prod = 30000  # kg
n_hours = 24

# 初始化全局模型
m = GEKKO(remote=False)
m.options.SOLVER = 1  # APOPT求解MINLP
m.solver_options = ['minlp_maximum_iterations 500',
                    'minlp_max_iter_with_int_sol 10',
                    'minlp_as_nlp 0',
                    'nlp_maximum_iterations 50',
                    'minlp_branch_method 1',
                    'minlp_integer_tol 0.05',
                    'minlp_gap_tol 0.01']

# 全局参数
c_H2 = m.Param(value=4.5, name="LCOH")
c_HB = m.Param(value=0.757, name="xCost_Haber-Bosch")
p_ue = m.Param(value=1.602, name="Price_europe")
p_usa = m.Param(value=1.323, name="Price_usa")
p_asia = m.Param(value=0.84, name="Price_asia")
p_br = m.Param(value=1.55, name="Price_brazil")
per_min_br = m.Param(value=20, name="minimmum_brasil")
perc = m.Param(value=100, name="total_percentage_ammonia")
limit = m.Param(value=963070.78)  # Electrolyzer limit (W)

# 时段参数
p_ufvi = m.Param(value=p_ufv.values, name="PV_power")
c_Ci = m.Param(value=c_C.values, name="grid_cost")

# 时段变量
per_eur = m.Array(m.Var, n_hours, lb=0, ub=100, name="ammonia_eur")
per_usa = m.Array(m.Var, n_hours, lb=0, ub=100, name="ammonia_usa")
per_asia = m.Array(m.Var, n_hours, lb=0, ub=100, name="ammonia_asia")
per_br = m.Array(m.Var, n_hours, lb=0, ub=100, name="ammonia_br")
p_C = m.Array(m.Var, n_hours, lb=0, ub=963070.78, name="grid_power")

# 逐时段约束
for t in range(n_hours):
    # 氨销售比例约束
    m.Equation(per_eur[t] + per_usa[t] + per_asia[t] + per_br[t] + per_min_br <= perc)
    # 电解器功率限制
    m.Equation(p_ufvi[t] + p_C[t] <= limit)

# 总生产约束：必须满足30吨配额
total_prod = m.Intermediate(m.sum([(p_ufvi[t] + p_C[t])/11085 for t in range(n_hours)]))
m.Equation(total_prod >= contract_prod)

# 总利润目标
total_profit = m.Intermediate(m.sum([
    (per_eur[t]*(p_ue - c_HB - c_H2/5.67) +
     per_usa[t]*(p_usa - c_HB - c_H2/5.67) +
     per_asia[t]*(p_asia - c_HB - c_H2/5.67) +
     (per_br[t]+per_min_br)*(p_br - c_HB - c_H2/5.67)) * (p_ufvi[t]+p_C[t])/1108500
    - c_Ci[t] * p_C[t]/1000000
    for t in range(n_hours)
]))

# 最大化总利润
m.Maximize(total_profit)

# 求解
try:
    m.solve(disp=True)
except Exception as e:
    print(f"求解失败: {e}")
    m.open_folder()

# 提取结果
nh3_ue = [v.value[0] for v in per_eur]
nh3_usa = [v.value[0] for v in per_usa]
nh3_asia = [v.value[0] for v in per_asia]
nh3_br = [v.value[0] for v in per_br]
p_buy = [v.value[0] for v in p_C]
prod_index = [(p_ufvi[t].value[0] + p_C[t].value[0])/11085 for t in range(n_hours)]
accumulate_prod = [sum(prod_index[:t+1]) for t in range(n_hours)]
max_profit = total_profit.value[0]

# 打印关键结果
print(f"总利润: {max_profit:.2f}")
print(f"总产量: {total_prod.value[0]:.2f} kg")
print(f"购电总量: {sum(p_buy):.2f} W")

3. 双目标的替代方案（可选） #

如果需要保留弹性产量要求，可以用目标规划给产量加一个远大于利润的权重，确保求解器优先满足产量：

# 目标规划：权重1e6远大于利润量级，优先满足产量
m.Maximize(1e6 * (total_prod - contract_prod) + total_profit)

但这种方式不如硬约束直观可靠，建议优先使用硬约束方案。

内容的