为何Xpress IVE线性规划模型忽略水力发电非负约束？

问题背景 #

我构建了一个以虚构企业能源产销利润最大化为目标的Xpress IVE线性规划模型，包含水力发电模块：

• 水库设有最大库容限制，存在自然降雨入流
• 低需求时段可消耗电力，以0.8的效率抽水入库
• 水库t时刻水量计算公式：t时刻水量 = t-1时刻水量 + 自然入流 + 抽水量 - 发电用水量

核心疑问 #

模型中已定义所有能源输出非负的约束：

NonNegativePowerCS(t, s) := PowerOutput(t, s) >= 0

但求解时却出现负水力发电值（模型以此实现水库蓄水），无法理解为何该约束被忽略。

模型代码 #

model EPower
  uses "mmxprs"
  
  declarations
  Emiss:                    set of string
    Times:                      set of real
    PeriodLength:             array(Times) of real
    Demand:                   array(Times) of real
    Source:                   set of string
    MaxOutput:                array(Source) of real
    RunningCost:              array(Source) of real
    IncreaseCost:             array(Source) of real
    Emissions:                array(Emiss, Source) of real
    EmissionsLimit:           array(Emiss) of real
    ElectricityPrice:         array(1..1) of real
    WaterInReserve:           array(Times) of mpvar ! Akin to store. This is ENERGY.
    WaterInStock:             array(Times) of mpvar ! Akin to stock. This is ENERGY.
    PowerToPumpUpWater:       array(Times) of mpvar ! Akin to make. This is POWER.
    MaxHydroReserve:          array(1..1) of real
    NaturalHydroInflow:       array(1..1) of real
    HydroPowerEfficiencyFrac: array(1..1) of real
    PowerOutput:              array(Times, Source) of mpvar
    PowerOutputChange:        array(Times, Source) of mpvar

end-declarations

initializations from "EPower.dat"
    Emiss Times PeriodLength Demand Source MaxOutput RunningCost IncreaseCost Emissions EmissionsLimit
    MaxHydroReserve NaturalHydroInflow HydroPowerEfficiencyFrac ElectricityPrice
end-initializations


! Power output constraints.
forall(t in Times, s in Source) do
  OutputCS(t, s) :=  PowerOutput(t,s) <= MaxOutput(s)
  NonNegativePowerCS(t, s) := PowerOutput(t, s) >= 0 ! WHY DOES THIS NOT WORK?????
end-do

forall(t in Times, s in Source) do
  if t = 1 then
    ChangeCS(t,s) := PowerOutput(t, s) = PowerOutputChange(t,s) + PowerOutput(6, s) ! These seem to work.
  else
    ChangeCS(t,s) := PowerOutput(t, s) = PowerOutputChange(t,s) + PowerOutput(t-1, s)
  end-if
end-do



forall(t in Times, s in Source) do
  ICPerTimePerSource(t, s) := PowerOutputChange(t, s) * IncreaseCost(s)
  ICPerTimePerSourceCS(t, s) := ICPerTimePerSource(t, s) >= 0 ! I think this is also a problem, but I'm not sure right now. Focussed on negative power.
end-do

forall(t in Times) do
  DemandCS(t) := sum(s in Source) PowerOutput(t, s) = Demand(t) + PowerToPumpUpWater(t) ! We need to generate enough power to meet demand and pump up water.
end-do

forall(t in Times, s in Source)
  EnergyProducedT(t, s) := PowerOutput(t, s) * PeriodLength(t) ! Total energy each source produces in each time period.

forall(s in Source)
  TotalEnergyProducedPerSource(s) := sum(t in Times) EnergyProducedT(t, s) ! Total energy each source produces.


! Emission constraints.
forall(s in Source, e in Emiss)
  EmissionsPerSource(e, s) := TotalEnergyProducedPerSource(s) * Emissions(e, s)

forall(e in Emiss) do
  TotalEmissions(e) := sum(s in Source) EmissionsPerSource(e, s)
  EmissionsCS(e) := TotalEmissions(e) <= EmissionsLimit(e)
end-do

! What follows are the hydro constraints.

forall(t in Times) do
    if t = 1 then
        WaterInStockCS(t) := WaterInStock(t) = WaterInStock(6) + PowerToPumpUpWater(t) * HydroPowerEfficiencyFrac(1) * PeriodLength(t) + NaturalHydroInflow(1) * PeriodLength(t) - PowerOutput(t, "Hydro") * PeriodLength(t) ! The problem is here.
        WaterInReserveCS(t) := WaterInReserve(t) = WaterInStock(6) + PowerToPumpUpWater(t) * HydroPowerEfficiencyFrac(1) * PeriodLength(t)  + NaturalHydroInflow(1) * PeriodLength(t)
    else
        WaterInStockCS(t) :=
          WaterInStock(t) = WaterInStock(t-1) + PowerToPumpUpWater(t) * HydroPowerEfficiencyFrac(1) * PeriodLength(t) + NaturalHydroInflow(1) * PeriodLength(t) - PowerOutput(t, "Hydro") * PeriodLength(t)
        WaterInReserveCS(t) :=
          WaterInReserve(t) = WaterInStock(t-1) + PowerToPumpUpWater(t) * HydroPowerEfficiencyFrac(1) * PeriodLength(t) + NaturalHydroInflow(1) * PeriodLength(t)
    end-if
    MaxReservoirCapacityCS(t) := WaterInReserve(t) <= MaxHydroReserve(1)
end-do

forall(t in Times) do
  WindCS(t) := PowerOutput(t, "Wind") <= 6000 ! This is only for stage 1.
end-do


TotalRunningCost := sum(s in Source) TotalEnergyProducedPerSource(s) * RunningCost(s)
TotalIncreaseCost := sum(t in Times, s in Source) ICPerTimePerSource(t,s)
TotalCost := TotalRunningCost + TotalIncreaseCost
TotalIncome := sum(s in Source) TotalEnergyProducedPerSource(s) * ElectricityPrice(1)
TotalProfit := TotalIncome - TotalCost

maximise(TotalProfit)

forall(s in Source) do
  !writeln("Water in reserve: ", getsol(WaterInReserve(t)))

  !writeln("Power output increase: ", getsol(PowerOutputChange(1, s)))
  writeln("Power output: ", getsol(PowerOutput(1, s)))
  !writeln("Water in stock: ", getsol(WaterInStock(t))
end-do
writeln("Water pump power: ", getsol(PowerToPumpUpWater(1)))
end-model

问题根源与修正方案 #

根源分析 #

1. 变量复用导致逻辑冲突
   你在水力平衡约束中，用- PowerOutput(t, "Hydro") * PeriodLength(t)表示发电耗水，但当PowerOutput为负时，这一项会变成增加水库水量——相当于把发电变量当成了蓄水变量使用，直接违背了PowerOutput的定义（发电输出），迫使模型违反非负约束来满足水力平衡。

2. 约束逻辑矛盾
   你已经定义了PowerToPumpUpWater(t)作为抽水用电变量，但未正确用它实现蓄水，反而试图通过反向发电替代，造成约束间逻辑矛盾。Xpress求解时会优先满足等式约束（水力平衡），而非不等式约束（非负输出）。

修正步骤 #

• 明确变量职责：PowerOutput(t, "Hydro")仅表示水力发电输出，必须严格非负；蓄水操作完全通过PowerToPumpUpWater(t)实现，对应的水量增加为PowerToPumpUpWater(t) * HydroPowerEfficiencyFrac(1) * PeriodLength(t)。
• 修正水力平衡约束：移除用负发电实现蓄水的逻辑，确保PowerOutput(t, "Hydro")在约束中仅作为耗水项（即减去正数值），蓄水仅通过PowerToPumpUpWater变量处理。
• 验证约束有效性：求解前添加display(NonNegativePowerCS)语句，检查所有非负约束是否被正确加载，避免因语法错误导致约束被忽略。

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