Why does my model ignore a linear constraint?

Question

I have an Xpress IVE linear programming model whose objective is to maximise the profit from producing and selling energy for a fictitious company. One of the sources of energy is hydro. There is a reservoir with max capacity, so within the model we also have an inventory problem, because there is a natural rain inflow into the reservoir, but at times of low demand we can also use power to pump up water to the reservoir at an efficiency of 0.8. So, roughly speaking, the water in the reservoir at time t is equal to the water in the reservoir at time t-1, plus the natural inflow, plus what we've pumped up, minus what we use to generate power.

The problem is this: as you shall see in the code below, I have a constraint that says specifically that the power output for any source must be non-negative. Yet the model is 'producing' negative hydro power because that means that it fills up the reservoir with it. So I'm not sure why my constraint is being ignored.

Here is the code:

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

Why does my model ignore a linear constraint?

Answers (1)

Related Questions